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Abstract 

High-Frequency  Acoustic  Volume  Scattering  from 
Biologically  Active  Marine  Sediments 

by  Christopher  D.  Jones 

Chair  of  Supervisory  Committee 

Professor  Darrell  R.  Jackson 
Electrical  Engineering 


A  thesis  on  high-frequency  acoustic  volume  scattering  from  marine  sediments 
with  application  to  remote  sensing  of  benthic  biological  activity  is  presented.  Small 
perturbation  theory  is  used  to  describe  bistatic  volume  scattering  in  a  sediment 
half-space.  The  sediment  is  modeled  as  an  acoustic  fluid  with  random  fluctuations 
in  density  and  compressibility.  Insight  into  determining  whether  single  or  multiple 
scattering  is  significant  in  the  medium  is  gained  by  using  the  bilocal  approximation 
to  Dyson’s  equation.  An  alternative  analysis  of  volume  scattering  is  made  using 
exact  numerical  simulations,  and  a  numerical  method  for  two-dimensional  volume 
scattering  using  the  method  of  moments  is  presented.  Both  periodic  and  nonperiodic 
random  media  are  considered.  Scattering  theory  is  compared  with  numerical  Monte- 
Carlo  simulations,  and  the  validity  of  the  small  perturbation  method  is  discussed. 
The  effects  of  the  half-space  scattering  geometry  on  the  coherent  field  within  the 
sediment  and  on  the  bistatic  scattering  cross-section  are  investigated. 

Benthic  biological  activity  creates  temporal  and  spatial  variations  in  the  sedi¬ 
ment  physical  properties  that  result  in  temporal  and  spatial  variations  in  sediment 


volume  scattering.  This  acoustic  variability  is  used  as  a  remote  sensing  tool  to  infer 
parameters  of  biological  activity,  or  bioturbation.  To  develop  a  forward  model  that 
relates  bioturbation  to  density  fluctuations  in  the  sediment  and,  therefore,  to  acous¬ 
tic  scattering,  a  new  stochastic  model  of  bioturbation  is  presented  that  describes 
biological  mixing  as  an  inhomogeneous  (two  scale)  biodiffusion  process.  Nonlocal 
mixing  (due  to  macrofauna)  is  described  as  a  filtered  Poisson  process,  and  local  mix¬ 
ing  (due  to  meiofauna)  is  described  as  diffusive.  Modeling  issues  such  as  the  spatial 
stationarity  of  bioturbation  are  discussed. 

The  bioturbation  and  acoustic  scattering  models  are  then  combined  to  produce 
a  model  for  the  decorrelation  in  time  of  acoustic  backscatter.  Model  predictions 
are  compared  with  experimental  data  collected  over  a  two  month  period  during 
the  Orcas  Island  experiment.  The  observed  decorrelation  of  acoustic  backscattering 
from  the  sediment  at  the  Orcas  site  is  compared  to  model  predictions  of  temporal 
decorrelation,  and  the  feasibility  of  using  acoustic  remote  sensing  to  detect  and  study 
benthic  biological  activity  is  discussed. 
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BENTHIC:  Pertaining  to  the  bottom  of  the  sea  or  lake. 

EPIFAUNA:  Benthic  organisms  that  live  on  the  surface  of  the  sediment-water 
interface. 
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MACROFAUNA:  Benthic  organisms  that  fall  with  the  size  category  of  approxi¬ 
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MEIOFAUNA:  Benthic  organisms  that  fall  within  the  size  category  of  approxi¬ 
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Chapter  1 

INTRODUCTION 

Acoustic  wave  propagation  and  scattering  have  been  widely  used  to  study  pro¬ 
cesses  and  detect  objects  in  the  ocean.  The  principal  reasons  are  the  desire  to 
monitor  and  collect  environmental  data  over  a  large  area  and  the  limited  usefulness 
of  optical  and  electromagnetic  wave  propagation  underwater.  For  a  wide  range  of 
scientific  and  engineering  problems  dealing  with  the  ocean  environment,  acoustic 
wave  propagation  and  scattering  is  the  only  remote  sensing  technology  available. 

The  ocean,  however,  is  a  highly  stochastic  environment.  For  many  practical 
applications  of  underwater  acoustics,  there  is  inadequate  oceanographic  information 
about  the  spatial  and  temporal  variability  in  the  medium  to  adequately  model  acous¬ 
tic  propagation  and  scattering.  In  some  cases,  even  the  mechanisms  by  which  sound 
is  scattered  are  not  well  understood  and  must  be  highly  idealized.  The  problem  is 
further  complicated  by  the  low  frequency  and  phase  velocity  of  acoustic  waves  (rela¬ 
tive  to  electromagnetic  waves).  In  many  cases,  the  ocean  environment  fluctuates  on 
spatial  scales  comparable  to  the  acoustic  wavelength,  especially  at  higher  frequen¬ 
cies.  Consequently,  the  modeling  of  underwater  acoustics  for  high-frequency  remote 
sensing  in  the  ocean  is  inherently  linked  to  stochastic  modeling  of  ocean  processes 
that  affect  propagation  and  scattering. 

In  general,  there  are  two  very  different  reasons  for  studying  acoustic  interaction 
with  the  ocean  environment.  One  reason  is  to  remove  the  negative  effects  of  the 
random  environment  on  acoustic  systems,  such  as  underwater  communication  de- 
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vices  and  object  detection  sonars.  The  goal  of  understanding  the  environment  is 
to  provide  better  forward  models  of  propagation  and  scattering  that  can  be  used  to 
improve  the  systems’  performance.  The  other  reason  is  to  facilitate  the  use  of  acous¬ 
tics  as  a  remote  sensing  tool  for  studying  the  ocean  environment.  A  comprehensive 
understanding  of  the  physics  of  scattering  in  the  ocean  will  lead  to  improved  inverse 
models  of  scattering  and  propagation.  With  reliable  acoustic  models,  over  a  range 
of  frequencies,  underwater  acoustics  can  be  used  to  study  the  ocean  on  spatial  and 
temporal  scales  that  are  unattainable  using  other  technologies. 

Some  aspects  of  the  ocean  acoustics  problem  have  been  investigated  for  many 
years,  especially  low  frequency  propagation  in  the  deep  water  environment.  In  re¬ 
cent  years,  considerable  attention  has  focused  on  coastal  environments  and  shallow 
water  processes,  and  the  use  of  high  frequency  acoustics  in  the  10-500  kHz  range. 
Propagation  and  scattering  at  these  frequencies  are  sensitive  to  small-scale  ocean 
processes  including:  air-sea  interactions  that  create  bubbles  and  surface  roughness; 
internal  wave  structures  and  turbulence  that  create  index  of  refraction  fluctuations 
in  the  water  column;  and  biological  and  hydrodynamic  processes  in  the  water  and 
sediment.  Seafloor  processes,  such  as  bioturbation  and  sediment  transport,  create 
sediment-water  interface  roughness  and  sediment  volume  heterogeneity. 

The  topic  of  this  thesis  is  the  use  of  high-frequency  acoustic  volume  scattering  in 
marine  sediments  as  a  remote  sensing  tool.  Two  issues  in  seafloor  scattering  will  be 
addressed  in  detail:  1)  the  more  general  issue  of  modeling  volume  scattering  from 
the  seafloor;  and  2)  the  specific  application  of  acoustic  remote  sensing  of  biological 
processes  in  the  seafloor  -  in  particular,  the  biological  processes  that  create  spatial 
and  temporal  variability  in  the  volume  of  the  sediment.  The  second  chapter  of 
this  thesis  will  address  the  problem  of  modeling  acoustic  volume  scattering  from 
a  fluid  half-space  random  medium.  Half-space  refers  to  the  geometry  created  by 

the  discrete  interface  between  the  sediment  and  the  overlaying  water.  The  small 

\ 

perturbation  method  will  be  used  to  formulate  analytic  expressions  for  the  first  and 
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second  moments  of  the  scattered  field  due  to  a  time-harmonic  (CW)  excitation. 

In  chapter  three,  a  numerical  method  for  volume  scattering  from  a  fluid  half-space 
will  be  presented.  The  numerical  solutions  are  exact  in  that  they  attempt  to  solve 
the  complete  integral  equation  formulation  of  the  problem,  including  attenuation  in 
the  medium,  multiple  scattering,  and  half-space  effects.  Approximations  to  the  wave 
solution  are  only  numerical.  With  increased  computational  resources,  the  numerical 
solution  will  improve  and  approach  the  exact  solution. 

Statistical  parameters  such  as  the  scattering  cross  section  and  field  correlation 
are  found  by  generating  random  realizations  of  sediment  heterogeneity  and  perform¬ 
ing  Monte-Carlo  simulations  of  the  scattered  field.  Only  the  time  harmonic  (CW) 
solution  is  considered.  Numerical  results  will  be  presented,  and  the  computational 
problems  and  convergence  issues  associated  with  performing  exact  two-dimensional 
calculations  of  volume  scattering  will  be  discussed. 

In  chapter  four,  the  validity  of  the  small  perturbation  method  for  scattering  from 
typical  sediments  will  be  discussed.  Multiple  scattering  theory,  in  the  context  of  the 
bilocal  approximation  to  Dyson’s  equation,  will  be  used  to  investigate  scattering 
due  to  density  and  compressibility  fluctuation.  Numerical  simulations  (guided  by 
the  multiple  scattering  theory)  will  also  be  used  to  investigate  multiple  scattering 
and  half-space  effects.  Numerical  results  will  be  compared  with  the  perturbation 
solution,  and  inferences  will  be  drawn  as  to  when  the  approximations  are  valid. 

The  last  two  chapters  (five  and  six)  comprise  the  second  issue  of  this  thesis:  an 
application  of  sediment  volume  scattering  to  the  remote  sensing  of  benthic  biolog¬ 
ical  activity.  In  chapter  five,  a  new  stochastic  model  of  bioturbation  is  presented 
in  which  the  spatial  and  temporal  fluctuations  of  sediment  physical  properties  are 
related  to  biological  activity.  Existing  models  of  sediment  diagenesis  are  reviewed 
and  the  mechanisms  by  which  biological  activity  creates  sediment  heterogeneity  are 
discussed.  Then,  in  chapter  six,  a  model  of  acoustic  scattering  due  to  bioturbation  is 
discussed.  By  coupling  the  sediment  bioturbation  model  with  the  scattering  model, 
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a  model  for  the  decorrelation  in  time  of  the  scattered  field  due  to  biological  activity 
is  developed.  This  model  is  then  compared  with  acoustic  field  data  and  sediment 
physical  data  obtained  in  a  biologically  active  shallow-water  environment. 


2*1  Volume  Scattering  in  Sediment 

An  acoustic  field  incident  upon  a  random  medium  will  interact  with  heterogeneities 
within  the  volume  to  create  a  scattered  field.  The  solution  for  the  volume  scat¬ 
tered  field  begins  by  examining  the  scales  and  interrelationships  between  various 
parameters  of  the  problem  at  hand.  The  most  essential  of  these  parameters  are: 
1)  the  magnitude  of  the  random  physical  properties  of  the  medium  that  give  rise 
to  scattering;  2)  the  characteristic  length  scales  of  the  randomness  in  the  medium 
relative  to  the  wavelength  of  the  incident  field;  and  3)  the  absorption  in  the  medium 
which  defines  the  characteristic  distance  over  which  scattering  occurs.  Typically, 
some  type  of  analytic  or  heuristic  criterion  is  applied  to  these  parameters,  and  a 
solution  strategy  is  chosen. 

There  are  two  general  approaches  to  the  solution  of  wave  scattering  in  a  random 
medium:  analytic  theory  and  radiative  transport  theory.  Analytic  wave  scattering 
theories  are  solution  methodologies  applied  directly  to  the  wave  equation  for  the 
field  quantity  of  interest.  They  can  be  broadly  categorized  as  either  single  or  mul¬ 
tiple  scattering  theories.  Radiative  transport  theory  is  based  on  phenomenological 
observations  of  the  transport  characteristics  of  the  wave  field  intensity  [17,  30].  The 
relationship  between  transport  theory  and  multiple  scattering  theory  has  been  well 
established  [80,  30].  However,  this  thesis  will  focus  mainly  on  analytical  approaches 
to  random  medium  scattering,  and  attention  will  be  given  to  the  criteria  for  de¬ 
termining  which  type  of  analytic  solution  methodology  is  applicable  to  sediment 
volume  scattering. 

A  single  (or  weak)  scattering  medium  is  one  in  which  the  total  perturbations 
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experienced  by  the  incident  field  are  small,  and  single  scattering  dominates  the  scat¬ 
tered  field.  Weak  scattering  methods  can  provide  analytical  solutions  to  problems 
where  general  solutions  are  not  attainable.  They  can  also  provide  intuition  into 
the  scattering  phenomena.  First-order  perturbation  methods,  such  as  the  Born  and 
Rytov  approximations,  are  classical  examples  of  weak  scattering  methods. 

Multiple  scattering  occurs  if  the  perturbation  of  the  incident  field  by  the  medium 
is  comparable  to  the  incident  field  itself  and  there  is  scattering  between  the  differ¬ 
ential  elements  of  the  volume.  Higher-order  perturbation  methods  can  be  used  to 
provide  insight  into  multiple  scattering.  In  general,  however,  they  are  complicated 
and  are  difficult  to  extend  to  higher  moment  solutions  of  the  scattered  field.  Effec¬ 
tive  medium  theories  can  be  used  to  define  renormalized  parameters  of  the  medium 
that  account  for  multiple  scattering.  For  example,  Dyson’s  equation  describes  the 
coherent  (mean)  field  within  a  scattering  medium.  It  can  be  used  to  define  a  com¬ 
plex  effective  wavenumber,  where  the  wave  speed  and  attenuation  include  the  effects 
of  multiple  scattering  [10,  75].  However,  effective  medium  methods  are  also  difficult 
to  extent  to  higher  moments  and  in  general  do  not  provide  analytic  solutions  for 
interesting  problems. 

Single  scattering  approximations  are  used  in  many  ocean  acoustics  applications. 
For  volume  scattering  from  marine  sediments,  the  Born  approximation  is  widely 
assumed  [34,  26,  88].  However,  the  validity  of  the  weak  scattering  assumption  for 
typical  marine  sediments  is  not  well  established.  Typically,  sediment  core  sampling 
[13,  34]  or  in  situ  analyses  [74,  87]  are  used  to  estimate  density  and  sound  speed 
fluctuations  in  the  sediment.  The  measured  sediment  fluctuations  are  then  used  to 
estimate  backscatter  strength  using  single  scattering  models. 

Experimentally  measured  backscatter  strengths  have  been  compared  with  the 
predictions  of  first-order  perturbation  models  for  volume  scattering  [34].  Some  data 
show  reasonable  agreement  with  models,  within  the  confidence  limits  of  the  esti¬ 
mated  sediment  parameters  used  as  inputs  to  the  models.  However,  the  results  are 


6 


inconclusive.  No  attempt  is  made  to  investigate  half-space  effects  or  higher-order 
scattering  effects.  The  disagreements  in  model-data  comparisons  are  attributed  to 
errors  and  resolution  problems  associated  with  the  core  analysis  and  the  spectral  es¬ 
timation  of  the  sediment  fluctuations.  Comprehensive  ground-truth  measurements 
of  sediment  volume  heterogeneity  in  a  wide  variety  of  environments  are  not  yet 
available. 

In  contrast,  scattering  from  the  sediment- water  and  air-water  interface  in  the 
ocean  has  been  studied  in  detail  and  the  parameter  space  of  validity  determined 
[76,  77].  It  is  well  established  that  small  perturbation  methods  for  sediment  interface 
roughness  scattering  are  not  always  applicable  for  the  range  of  acoustic  frequencies, 
grazing  angles,  and  interface  statistics  encountered  in  the  ocean.  In  general,  a  variety 
of  surface  scattering  models  (e.g.,  perturbation,  tangent-plane,  small-slope  methods) 
are  applied  depending  on  the  particular  merits  of  the  method  and  the  characteristics 
of  the  surface  roughness  [37,  15].  Much  of  this  evaluation  was  done  in  comparison 
with  exact  numerical  solutions,  the  most  common  calculation  being  the  computation 
of  the  scattering  integral  equation  of  an  idealized  one-dimensional  rough  surface  by 
quadrature  methods  [76,  77].  More  recently,  efficient  numerical  techniques  have  been 
extended  to  the  solution  of  the  full  three-dimensional  surface  problem,  [78]  and  [22] 
for  example. 

This  thesis  will  address  the  validity  of  small  perturbation  methods  for  volume 
scattering  in  marine  sediments  (where  the  magnitudes  of  the  medium  fluctuations 
are  the  small  parameters).  Insight  into  determining  whether  single  or  multiple  scat¬ 
tering  is  significant  in  the  sediment  will  be  gained  by  using  classical  multiple  scat¬ 
tering  theory.  Effective  medium  theory  using  the  bilocal  approximation  to  Dyson’s 
equation  will  be  extended  to  include  density  fluctuation  (analogous  with  including 
permeability  fluctuation  in  the  EM  case).  An  alternative  analysis  (guided  by  the 
results  of  the  multiple  scattering  theory)  will  be  made  using  “exact”  numerical  sim¬ 
ulations  of  the  scattered  field.  By  comparing  solutions  of  the  perturbation  method 
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with  the  numerical  solutions,  the  range  of  validity  of  the  perturbation  methods  for 
typical  sediments  can  be  inferred.  Parameters  of  the  sediment  are  limited  to  those  of 
typical  marine  sands  and  muds  modeled  as  acoustic  fluid  half-spaces  with  Gaussian 
random  fluctuations  in  density  and  compressibility. 

Numerical  simulations  will  be  limited  to  the  solution  of  the  volume  integral 
equation  in  two  dimensions  and  for  a  single  frequency  using  the  method  of  mo¬ 
ments  [24].  Moment  methods  are  commonly  applied  to  electromagnetic  scattering 
problems  (see  [18]  and  [57]  for  an  overview).  Other  solution  techniques  such  as 
finite- volume  methods,  finite-difference  and  finite  element  methods  are  also  applied 
to  wave  propagation  and  scattering  (also  see  [57]  for  an  overview  of  these  and  other 
methods).  The  numerical  integral  equation  formulation  provides  one  obvious  ad¬ 
vantage:  it  is  the  most  direct  and  analogous  solution  methodology  for  comparison 
with  the  perturbation  results.  By  treating  the  scattering  integral  as  a  matrix  equa¬ 
tion,  the  convergence  properties  of  the  perturbation  expansion  can  be  studied  by 
examining  the  convergence  properties  of  the  solution  of  the  matrix  equation. 


1.2  Benthic  Biological  Activity 

Biologically  active  sediments  are  continually  modified  by  epifauna  and  infauna.  Epi- 
fauna  are  organisms  that  live  on  the  sediment  surface,  and  infauna  are  organisms 
living  within  the  sediment.  The  collective  mixing  effect  of  these  organisms  is  termed 
bioturbation.  Interests  in  bioturbation  range  from  applications  in  marine  engineering 
to  biological  oceanography  including  the  time  evolution  of  acoustic  and  structural 
properties  of  seafloor  sediments,  and  the  monitoring  of  the  effects  of  pollutants  or 
other  disturbances  on  benthic  ecology.  Biological  interests  include  diagenesis  in  ma¬ 
rine  sediments,  the  study  of  sediment  disturbance  processes,  deposit-feeding  levels, 
and  population  dynamics,  for  example. 

Observations  of  the  influence  of  benthic  biology  on  sediment  heterogeneity  are 
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Figure  1.1:  Examples  of  benthic  epifauna  and  infauna  that  create  bioturbation  in 
marine  sediments  (from  Berner  1980) 

numerous.  An  overview  of  this  topic  can  be  found  in  texts  on  the  subject  [4,  40,  46] 
and  in  review  articles  with  extensive  bibliographies  [56],  for  example.  Figure  1.1 
illustrates  typical  organisms  of  interest  in  this  development.  Sediment  x-ray  cores 
reveal  the  structure  of  density  fluctuations  (through  x-ray  absorption)  within  the 
volume  of  the  sediment.  Figure  1.2  shows  several  examples  of  typical  X-radiographs 
taken  in  biologically  active  sediments.  The  imagery  shows  typical  burrows  and 
inclusions  due  to  biological  activity.  These  illustrations  and  images  also  clarify  the 
spatial  scales  of  scattering  that  are  of  primary  interest  in  this  thesis.  The  acoustic 
wavelength  for  frequencies  of  interest  (10-500  kHz  @@  1500  m/s)  is  comparable  to 
the  scale  of  the  heterogeneity  created  by  biological  processes. 

Acoustic  observations  of  benthic  biological  activity  are  also  numerous,  although 
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the  observations  are  not  generally  investigated  for  the  purpose  of  understanding 
the  biota.  Of  particular  interest  are  recent  experiments  in  which  observations  of 
sediment  backscatter  were  made  over  a  period  of  months  in  a  variety  of  shallow- 
water  environments  and  sediment  types  [20,  36,  41,  86].  The  spatial  and  temporal 
evolution  of  scattering  from  the  seafloor  (at  40  and  300  kHz)  has  been  attributed  to 
biological  reworking  of  sediment  microtopography,  sediment  volume  heterogeneity, 
and  near-bottom  volume  reverberation.  The  correlation  in  time  of  the  backscattered 
field  shows  characteristic  decay  that  is  most  likely  a  function  of  the  type  and  rate 
of  biological  activity  associated  with  the  different  environments  (see  Figure  1.3). 
During  the  Orcas  experiment,  changes  in  the  acoustic  signatures  of  buried  objects 
over  time  were  observed  [86].  Biological  activity  around  the  buried  objects  decreased 
the  object’s  total  scattering  strength  in  time.  During  the  same  time  period,  the 
temporal  correlation  of  backscatter  from  the  area  of  the  object  showed  increased 
decorrelation  -  most  likely  due  to  increased  biological  activity  around  the  object. 

On  micro-  and  macroscopic  scales,  bioturbation  affects  the  physical  and  chemical 
properties  of  the  sediment.  The  intensity  of  bioturbation  affects  the  distribution  of 
marine  chemicals  which,  in  turn,  influence  the  microbial  activity*  in  the  sediment 
and  the  ultimate  fate  of  marine  pollutants  [1,  4,  40,  83].  Epifaunal  activity  creates 
microtopography  that  may  decrease  the  critical  velocity  necessary  to  erode  surface 
layers  [84].  Infaunal  activity,  such  as  tube  building,  is  responsible  for  vertical  and 
horizontal  redistribution  of  solid  material  within  the  sediment,  creating  spatial  and 
temporal  heterogeneity  in  sediment  porosity  [6,  7]  which  leads  to  heterogeneity  in 
the  sediment  bulk  density  and  compressibility.  Often,  the  traces  of  these  mixing 
activities  provide  the  only  evidence  of  the  existence  of  fauna  in  an  area  and  may 
provide  a  measure  of  benthic  biomass  [84] . 

Bioturbation  is  pervasive  and  shows  similar  spatial  characteristics  across  a  wide 
variety  of  marine  environments.  The  effective  biological  mixing  depth  (based  on  a 
variety  of  measurement  methods)  has  a  worldwide  mean  of  approximately  10  cm 
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Figure  1.3:  Decorrelation  of  the  backscattered  field  due  to  biological  activity  at 
various  shallow  water  locations. 


12 


Figure  1.4:  Spectrum  of  the  equivalent  spherical  diameter  body  size  of  benthic 
organisms  in  marine  sediments  (from  Jumars  1993). 


[40].  In  addition,  the  body  size  spectra  of  benthic  organisms  show  characteristic 
tri-modal  structure,  rather  than  a  flat  spectrum  observed  for  organisms  in  the  ocean 
water  column  [40].  Figure  1.4  shows  the  typical  size-abundance  relationship  for 
organisms  in  marine  sediments.  In  both  shallow  water  and  deep  sea  sediments, 
there  are  minima  in  the  body  size  spectra  at  about  10  //m  and  1  mm.  In  general, 
meiofauna  collectively  outweigh  macrofauna,  with  animal  and  bacterial  abundance 
declining  exponentially  as  a  function  of  depth  into  the  sediment  [40] . 

This  brief  introduction  is  given  mainly  to  motivate  interest  in  bioturbation  from 
a  remote  sensing  point  of  view,  not  as  a  comprehensive  review  of  the  topic.  The 
subject  of  bioturbation  in  marine  sediments  is  far  beyond  the  scope  and  intent  of 
this  thesis.  However,  because  one  of  the  objects  for  this  thesis  is  to  address  the 
remote  sensing  problem,  a  forward  model  of  bioturbational  mixing  must  be  devel¬ 
oped.  The  forward  model  in  this  case  is  a  quantitative  description  of  mixing  that 
can  relate  biological  parameters  to  acoustically  detectable  changes  in  the  physical 
properties  of  the  sediment.  A  certain  level  of  detailed  description  of  the  biological 
processes  is  necessary,  therefore,  to  justify  the  proposed  modeling  methodology.  Fur¬ 
thermore,  since  the  mixing  model  is  used  as  the  critical  link  between  the  observation 
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of  sound  scattered  from  the  sediment  and  the  estimation  of  parameters  that  describe 
the  sediment  biology,  it  should  incorporate  established  variables  and  concepts  used 
by  benthic  biologists.  In  addition,  the  mixing  model  should  provide  the  physical 
reasoning  that  relates  these  parameters  to  the  modeling  of  sound  scattering. 
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Chapter  2 

SEDIMENT  VOLUME  SCATTERING 

In  general,  acoustic  scattering  from  the  seafloor  is  a  function  of  the  roughness 
at  the  sediment-water  interface,  roughness  of  interstitial  sediment  layers,  the  distri¬ 
bution  of  discrete  particles  or  bubbles  within  the  sediment,  and  fluctuations  in  the 
sediment  volume  density  and  sound  speed.  The  sediment  itself  can  be  modeled  as 
a  fluid  [32,  34],  an  elastic  [35]  or  poro-elastic  medium  [72,  19],  or  as  a  system  of 
densely  packed  sediment  grain  particles.  The  variety  in  model  idealizations  is  due  in 
part  to  the  wide  variety  of  sediment  types  and  seafloor  compositions  encountered  in 
the  ocean  environment  and  the  wide  range  of  acoustic  frequencies  used  to  interact 
with  the  seafloor. 

In  any  particular  scenario  of  sound  interaction  with  the  seafloor,  one  or  more 
of  the  above  mentioned  scattering  mechanisms  or  model  idealizations  may  be  ap¬ 
plicable.  Typically,  one  scattering  mechanism  is  assumed  to  be  dominant.  Surface 
roughness  scattering  is  likely  the  dominant  scattering  mechanism  for  many  hard 
sand  and  rock  sediments.  Volume  scattering  is  dominant  in  sediments  composed 
of  soft  silts  and  mud  [33].  In  the  case  of  competing  scattering  mechanisms,  where 
two  or  more  types  of  scattering  are  present  and  of  comparable  effect,  single  scatter¬ 
ing  is  typically  assumed  and  each  scattering  mechanism  treated  independently.  This 
chapter  will  address  the  issues  of  modeling  volume  scattering  from  marine  sediments. 

The  appropriate  medium  idealization  to  apply  at  a  particular  frequency  (inde¬ 
pendent  of  the  scattering  mechanism)  is  still  a  matter  of  some  debate.  In  this  thesis, 
it  will  be  assumed  that  sediment  can  be  modeled  as  an  effective  fluid  medium.  No 
statement  is  made  about  the  sediment  particle-wave  interaction  at  the  frequencies 
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Figure  2.1:  Geometry  for  volume  scattering  from  a  half-space  random  medium. 


of  interest,  only  that  the  net  result  is  an  effective  continuum  whose  acoustic  proper¬ 
ties  are  modeled  as  a  fluid.  Scattering  from  larger  particles  such  as  shells  and  drop 
stones  will  not  be  considered.  It  is  further  assumed  that  the  compressional  wave 
attenuation  of  the  sediment  is  high  and  acoustic  penetration  is  limited  to  the  region 
of  sediment  near  the  sediment-water  interface.  For  highly  consolidated  sediments, 
shear  properties  may  be  significant  and  the  medium  would  be  modeled  as  an  elastic 
solid.  However,  in  soft  sand  or  mud  (near  the  sediment-water  interface)  the  acoustic 
shear  properties  of  the  medium  are  weak,  and  the  sediment  will  be  treated  as  an 
acoustic  fluid  [35]. 

Figure  2.1  illustrates  the  fluid-fluid  half-space  interface  and  geometry  for  volume 
scattering.  The  grazing  and  scattering  angles  are  defined  relative  to  the  plane  of 
the  interface.  The  upper  half-space  (UHS)  contains  the  positions  of  the  source  and 
receiver  and  is  assumed  to  be  acoustically  homogeneous.  The  lower  half-space  (LHS) 
is  a  continuous  random  medium  of  infinite  extent  in  depth  and  in  the  horizontal 
dimension. 
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Figure  2.2:  Composition  of  marine  sediment  near  sediment-water  interface. 


2.1  Sediment  Physical  Properties 

The  sediment  is  assumed  to  be  unconsolidated  and  composed  of  a  mixture  of  solid 
sedimentary  particles  (or  grains)  and  sea  water  Figure  2.2  illustrates  a  typical  sed¬ 
iment  composition  near  the  sediment-water  interface.  Sediment  grains  may  consist 
of  a  mixture  of  silicate  sand  (quartz),  shell  fragments  (CaCOs),  clay  particles,  and 
organic  matter,  for  example.  The  sediment  is  assumed  to  possess  no  rigid  frame  (the 
shear  modulus  of  the  medium  is  zero)  and  is  treated  as  a  fluid  half-space  continuum. 
The  acoustic  phase  velocity  (c)  is  defined  in  terms  of  the  sediment  bulk  density  ( p ) 
and  compressibility  (k): 

(2.1) 

Randomness  in  the  sediment  is  characterized  by  fluctuations  in  the  sediment 
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porosity,  which  lead  to  fluctuations  in  the  sediment  density  and  compressibility. 
Normalized  fluctuations  in  the  sediment  density  and  compressibility  are  defined  as 

7 '«  =  («-  k)/k  , 

7 p  =  (P~  p)/p  > 

where  p  and  k  are  the  mean  values.  In  general,  the  sediment  density  is  a  real  quan¬ 
tity,  and  the  compressibility  is  complex  representing  dissipation  of  acoustic  energy. 
They  combine  to  create  a  complex  compressional  phase  speed  in  the  sediment.  The 
fluctuations  in  the  compressibility  are  assumed  to  have  the  same  relative  loss  as  the 
mean,  or 

Im  [£]  =  0  ;  (2.3) 

therefore  /yK  is  always  real.  The  fluctuation  in  sediment  sound  speed  is  defined  as 

7 c  =  (c-  c)/c  ,  (2.4) 

where 

(2.5) 

The  sound  speed  fluctuations  can  be  written  in  terms  of  the  density  and  compress¬ 
ibility  fluctuations: 


In  the  water,  p0  and  k0  are  real.  The  mean  values  of  density  or  compressibility  in 
both  the  UHS  and  LHS  are  assumed  to  be  constant  in  space  and  time  (no  gradients). 

In  general,  the  real  and  imaginary  components  of  the  wavenumber  may  take  any 
functional  form  that  satisfies  causality  (i.e.,  they  obey  the  Kramers-Kronig  relation¬ 
ship).  For  convenience,  however,  the  components  are  taken  to  be  proportional,  such 
that  the  wavenumber  in  the  LHS  is  defined  as 


k  =  Re[u;/c](l  +  iS )  , 


(2.7) 
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where  u  is  the  acoustic  frequency.  The  compressional  loss  factor  S  is  the  ratio  of 
the  imaginary  to  real  components  of  the  LHS  wavenumber.  It  is  assumed  to  be 
independent  of  frequency.  This  assumption  violates  causality  in  the  strict  sense. 
However,  over  the  frequency  range  of  interest  the  resulting  dispersion  errors  are 
small  and  can  be  neglected  [47]. 

If  S  is  independent  of  frequency,  and  dispersion  in  the  medium  is  negligible,  the 
attenuation  will  scale  as  the  first  power  of  frequency.  Extensive  experimental  evi¬ 
dence  suggests  the  first  power  frequency  dependence  [23,  60],  although  the  frequency 
dependence  of  attenuation  in  unconsolidated  marine  sediment  is  still  a  matter  of  de¬ 
bate.  Intergranular  friction  and  hysteresis,  for  example,  have  been  proposed  as  a  loss 
mechanism  in  which  linear  scaling  with  frequency  and  weak  logarithmic  dispersion 
can  be  explained  [16].  Porous  medium  (Biot)  models  with  added  dissipation  terms 
have  also  been  proposed  to  explain  approximately  linear  frequency  relations  [72,  19]. 

The  attenuation  coefficient  in  the  sediment  due  to  intrinsic  absorption,  with  units 
of  decibels  per  unit  wavelength  (dB/A),  is  defined  as 

407T  „ 


OL\  = 


ln(10) 


(2.8) 


Another  attenuation  coefficient  that  is  of  practical  use  has  units  of  decibels  per  meter 
per  kHz  (dB/m/kHz),  and  is  defined  as 


1000 

£*/  =  —Oi\  . 

c 


(2.9) 


A  constitutive  relationship  between  density  and  compressibility  is  not  available 
for  marine  sediments  composed  of  particles.  However,  a  linear  relation  is  assumed 
based  on  measured  density-sound  speed  relationships  [13,  45,  88,  32]: 


7*  =  VIp  >  (2-10) 

where  fi  is  a  constant,  at  least  within  the  confidence  intervals  of  the  available  data. 
Fluctuations  in  sediment  sound  speed  are  generally  weak  compared  to  fluctuations 
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in  density  and  compressibility,  indicating  a  negative  correlation  between  density  and 
compressibility.  A  negative  fj,  also  corresponds  to  a  physical  picture  of  sediments 
where  porosity  dictates  the  density  and  compressibility.  As  porosity  increases,  the 
bulk  density  of  the  sediment  decreases  (assuming  the  sediment  grains  are  more  dense 
than  water).  As  the  volume  of  water  increases  relative  to  the  volume  of  grains,  the 
compressibility  increases  (assuming  that  water  is  more  compressible  than  the  grains). 
As  density  increases,  compressibility  decreases.  Therefore,  fj,  is  generally  assumed 
to  be  negative. 

2.1.1  Correlation  and  Scale  of  the  Medium  Fluctuations 

Scattering  of  sound  from  the  sediment  volume  depends  on  the  spatial  correlation  of 
the  density  and  compressibility  fluctuations.  The  auto-  and  cross-correlations  of  the 
fluctuations  are  defined  as 

Cw,(ri,r2)  =  (7«(rx)7«(r2))  , 

CKp(rur2)  =  (7/t(ri)7/,(r2))  ,  (2.11) 

CV,(ri,r2)  =  (7P(rib,(r2))  • 

For  a  wide-sense  stationary  process,  the  correlation  is  defined  in  terms  of  the  differ¬ 
ence  coordinate  =  rx  —  r2: 

CKK(rd)  =  (7/e(r  +  rd)7/t(r))  .  (2.12) 

Similarly  for  the  other  correlation  functions.  If  it  is  assumed  that  the  medium 
compressibility  is  proportional  to  the  density  using  (2.10),  then 

CKp{rd)  =  CpK(rd )  =  fiCpp(rd)  ,  (2.13) 

and 


C««( rj)  =  ii2Cpr(rd)  . 


(2.14) 
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Various  functional  forms  of  the  correlation  function  can  be  assumed  depending 
on  the  spatial  scales  of  interest.  For  the  acoustic  frequencies  of  interest  and  the 
corresponding  scales  of  the  randomness  in  the  sediment  that  interact  with  the  high 
frequency  acoustic  waves  (correlation  scales  of  approximately  0.01-1.0  meters),  the 
von  Karman  or  the  exponential  correlation  functions  are  often  used  [30,  2,  69].  On 
scales  much  smaller  than  the  wavelength,  the  sediment  is  better  described  as  a  dense 
distribution  of  grain  particles,  and  a  description  in  terms  of  a  continuous  correlation 
function  is  not  applicable.  On  the  other  end  of  the  scale,  acoustic  interaction  with 
geological  features  (very  large  spatial  frequencies)  can  be  considered  deterministic 
and  handled  in  the  mean  sediment  properties.  For  example,  large-scale  sediment 
stratigraphy  is  observed  on  the  scale  of  meters  in  deep  cores.  However,  a  high- 
frequency  acoustic  field  will  only  penetrate  a  short  distance  into  the  sediment  due 
to  attenuation,  and  it  will  only  interact  with  a  single  layer. 

For  an  isotropic  medium,  the  von  Karman  correlation  function  is  defined  in  terms 
of  the  variance  of  the  density  fluctuations  (cr2p)  and  the  correlation  length  scale  of 
the  fluctuations  (/c)  as 

=  v(q-l)\K‘‘{rllc)  ’  (2'15) 

where  r  =  |r^|,  and  Kq  is  the  modified  Bessel  function  of  the  second  kind  of  order  q. 
The  exponential  correlation  function  is  a  special  case  of  the  von  Karman  correlation 
function  with  q  =  2: 

Cpp{rd)  =  cr2pe~r/lc  .  (2.16) 

The  exponential  correlation  function  is  a  commonly  assumed  form  for  sediment  core 
analysis  [34,  13].  The  main  reason  for  applying  the  exponential  assumption  is  more 
a  function  of  necessity  when  processing  the  sparsely  sampled  cores  rather  than  a 
statement  of  an  underlying  physical  process  in  the  sediment. 

The  power  spectral  density  (PSD)  function  of  a  stationary  process  is  defined  as 
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the  Fourier  transform  of  the  correlation: 

Sep(k)  =  C„(r)e-ik"*  .  (2.17) 

The  corresponding  power  law  PSD  function  for  the  exponential  correlation  function 
is 


s  (1-)- 

pA  )  (1//2  +  At2)2  * 

For  a  two  dimensional  process,  the  power  law  PSD  function  is 

PpK  )  (l//c2  +  fc2)3/2  • 


(2.18) 


(2.19) 


2.2  Volume  Scattering  Integral  Equation 


A  pressure  wave,  incident  upon  the  LHS,  will  interact  with  heterogeneity  in  the 
volume  to  create  a  scattered  field.  In  the  absence  of  fluctuations  in  the  medium,  the 
time-harmonic  wave  equation  valid  in  the  UHS  is 


V2p(r)  +  fcoP(r)  =  0  , 


(2.20) 


and  valid  in  the  LHS  is 


Pi r)V  • 


1 

Ar) 


Vp(r) 


+  k2p(r)  =  0  . 


(2.21) 


The  boundary  conditions  at  the  interface  between  the  two  fluids  (at  z  =  0)  are  the 
combined  Dirichlet  and  Neumman  conditions  for  a  penetrable  medium, 


P(x,y,z  =  0+)  =p(x,y,z  =  0  ) 


(2.22) 


and 


1  dp[x,  y,  2  =  0+)  _  1  dp(a:,  y,  £  =  0  ) 
Po  dz  p  dz 


(2.23) 
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In  the  presence  of  heterogeneity  in  the  LHS,  the  coefficients  of  the  wave  equation 
are  random  variables, 

/>(r)  V •  f 4>VP(r)l  +  k2(r)p(r)  =  0  ,  (2.24) 

LMr)  J 

where  k2(r)  =  cu2p(r)/c(r)  is  the  random  LHS  wavenumber.  A  single  wave  equation, 
valid  in  the  LHS  (with  or  without  heterogeneities),  can  be  defined  by  rearranging 
(2.24)  and  grouping  the  random  terms  on  the  right  side: 

V2p(r)  +  Pp( r)  =  -/( r)  .  (2.25) 

The  right  side  represents  an  induced  source  term,  /,  which  combines  scattering  due 
to  fluctuations  in  the  medium’s  compressibility  and  density: 

f(r)  =  Pn(r )p(r)  -  V  •  [7,(r)Vp(r)]  .  (2.26) 

The  left  side  of  (2.25)  is  the  unperturbed  Helmholtz  operator.  Equation  (2.25) 
is  still  homogeneous  in  the  unknown  quantity  p.  No  new  energy  is  added  into  the 
system  by  the  source  term.  The  sources  act  to  redistribute  the  incident  wave  through 
interactions  with  the  heterogeneity  in  the  medium. 

The  Green’s  function  for  the  unperturbed  system  is  the  solution  of  the  wave 
equations 

V3s»(r)  +  P*,(r)  =  S(r  -  r')  (2.27) 

in  the  LHS,  and 

V2<7o(r)  +  fco<7o(r)  =  — £(r  —  r')  (2.28) 

in  the  UHS,  that  satisfies  the  appropriate  boundary  conditions  along  the  sediment- 
water  interface  (see  Appendix  A).  With  no  gradients  in  the  mean  density,  the 
spectral  representation  of  the  half-space  Green’s  function  for  transmission  from  r' 
in  the  LHS  to  r  across  the  interface  into  the  UHS  is 

OO 

<7o(r,r')  =  - ^  J  > 


(2.29) 
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where 


p  =  Jw-ki-ki 


(2.30) 


h  =  A2  -  h  - 


(2.31) 


are  positive  to  represent  a  physical  field.  The  coefficient  T'  is  the  plane-wave  pressure 
transmission  coefficient  from  the  LHS  to  the  UHS: 

T'  =  lPof3_0  .  (2.32) 

Pofi  +  pPo 

For  r  and  r'  both  in  the  LHS,  the  half-space  Green’s  function  including  reflection  is 

OO 

go( r,r')=g^  J  ^  [R(kx,ky)ei0]z+z'\  +  e^z-z'^  eik^x~x'^ik^y-y,Ukxdky  ,  (2.33) 


where  R  is  the  plane- wave  pressure  reflection  coefficient: 

R  =  pofl  -  .  (2.34) 

pop  +  pft 0 

The  volume  scattering  integral  equation  is  obtained  by  applying  the  Kirchhoff- 
Helmholtz  integral  relationship  with  the  source  term  (2.26): 


0p(/)  dg0(r,r')l 

j  dn'  P[rJ  Bn'  dS 


+  f  f(r')go{r,r')dr' . 
Jv 


(2.35) 


The  first  integral  in  (2.35)  is  a  surface  integral  over  a  bounding  sphere  at  infinity, 
and  simply  reduces  to  the  incident  field.  The  second  integral  is  the  superposition 
and  propagation  of  the  induced  sources  with  the  unperturbed  Green’s  function.  The 
volume  integral  can  be  simplified  by  integrating  by  parts  and  applying  the  divergence 


theorem  such  that 


Jv  V'  •  hv(r')V'Kr')] go(r,  r')dr'  =  jf  jp(r')?^-go(r,  r *)ds' 

-  [  7p(r,)V/go(r,  r')  •  V'p(r')dr'  . 

Jv 


(2.36) 
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By  applying  the  radiation  condition  at  infinity,  the  surface  integral  in  (2.36)  vanishes. 
The  final  integral  equation  for  the  total  field  can  be  written  as  the  sum  of  the 
unperturbed  field  (po)  and  the  scattered  field  (ps): 

p(r)  =  pQ(r)  +  ps(r)  .  (2.37) 

The  scattered  field  is  the  sum  of  all  the  fields  due  to  the  differential  sources  within 
the  volume: 

ps( r)  =  /  [k2iK(r')g0(r,  r')p(r')  -f  7(5(r')V,£f0(r,  r')  •  V'p(r')]  dr'  .  (2.38) 

Jv 

Equation  (2.38)  is  similar  to  the  well-known  integral  equation  [54]  which  uses 
the  free-space  Green’s  function.  However,  now  the  half-space  boundary  conditions 
are  included  by  the  use  of  the  half-space  Green’s  functions.  The  first  term  in  the 
volume  integral  represents  monopole  scattering  due  to  compressibility  fluctuations 
at  a  point  in  the  volume.  The  second  term  represents  dipole  scattering  due  to  the 
dot  product  of  gradients  in  the  field  and  Green’s  function  at  a  point  in  the  volume. 
The  integration  is  over  the  volume  of  the  LHS  medium.  With  significant  attenuation 
in  the  lower  medium,  the  scattered  field  will  penetrate  only  a  limited  distance  into 

t 

the  medium  so  that  the  acoustic  energy  from  the  incident  field  will  be  confined  to  a 
layer  near  the  surface  of  the  LHS.  Therefore,  the  depth  dimension  of  the  volume  is 
finite  and  for  many  practical  problems  will  be  relatively  small. 

If  the  medium  is  fluctuating  in  time,  the  scattered  field  will  fluctuate  in  time, 
as  well  as  in  space.  If  the  medium  is  fluctuating  slowly  compared  to  the  frequency 
of  the  acoustic  wave,  the  time  dependence  of  the  medium  carries  through  into  the 
scattering  integral  equation: 

ps(r,t)=  [  [p7«(r,,%o(r,r>(r,,f)  + 

Jv  (2.39) 

7,(r',*)V'$b(r,r')  •  S7’p(r',t)]dr'  . 

Here,  the  time  dependence  is  a  function  of  the  medium  only.  The  pressure  field  is 
still  assumed  to  be  time-harmonic,  and  this  time  dependence  is  suppressed. 
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If  the  medium  density  and  compressibility  are  known,  the  scattered  field  inside 
and  outside  the  scattering  medium  can  be  found  uniquely  using  the  appropriate 
Green’s  function.  However,  analytic  solutions  to  (2.38)  and  (2.39)  for  a  random 
medium  are  generally  not  available. 

2.3  Small  Perturbation  Method 

If  the  fluctuations  in  the  medium  are  small,  the  solution  for  the  perturbed  system 
will  be  similar  to  the  solution  of  the  unperturbed  system.  Therefore,  assume  that 
the  unknown  scattered  field  can  be  expanded  as  a  perturbation  series  in  powers  of 
some  small  and  dimensionless  quantity  e,  such  that 

p(r)  =  (j)0  +  (r)  +  e202(r)  +  ■  •  ■  ,  (2.40) 

where  e  is  of  the  order  of  magnitude  of  the  medium  fluctuations: 

0(e)  =  0( 7«)  =  0(7P)  .  (2.41) 

The  perturbation  series  is  defined  as  an  asymptotic  series  with  respect  to  the  se¬ 
quence  1,  e,  e2,  •  •  •  and  the  basis  functions  4>n  having  the  property 

N 

P~Yj  =  °(eN)  as  e  -+  0  .  (2.42) 

n~0 

This  property  does  not  guarantee  that  the  series  will  converge.  Therefore,  increasing 
the  number  of  terms  in  (2.40)  does  not  necessarily  improve  the  solution.  Even  if  the 
series  does  seem  to  converge,  it  may  not  necessarily  converge  towards  the  correct 
result.  The  reason  for  this  is  that  the  asymptotic  expansion  only  makes  a  statement 
about  the  series  in  the  limit  as  e  ->  0,  not  as  N  — ¥  oo.  Therefore,  as  the  small 
parameter  e  increases  (as  the  medium  fluctuations  grow  larger),  the  perturbation 
solution  is  unpredictable.  However,  if  the  parameter  e  is  “small”,  the  higher-order 
terms  in  the  series  will  often  decay  rapidly  [42].  Therefore,  a  relatively  small  number 
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of  terms  can  be  used  to  accurately  represent  the  behavior  of  the  solution.  Each 
higher-order  term  in  the  series  then  corresponds  to  a  higher-order  correction  to  the 
true  field. 

The  consecutive  terms  of  the  perturbation  series  are  obtained  by  substituting 
(2.40)  into  the  integral  equation  (2.37)  and  balancing  terms  by  powers  of  e: 

P  =  fa  +  e<f>i  +  e2<f> 2  4 -  (2  43) 

=  po  4*  L[fa]  +  L[e(f)  x]  +  L[e2  fa]  +  •  •  •  . 

The  operator  L  is  the  linear  integral  operator  of  (2.38)  or  (2.39).  Since  only  single 
powers  of  7*  or  7P  are  present  in  the  operator,  L  is  first-order  in  e.  The  zeroth-order 
term  (corresponding  to  no  randomness  in  the  medium)  is  found  to  be 

fa  =  Po  •  (2.44) 

This  is  the  known  solution  to  the  unperturbed  problem,  and  for  a  unit  point  source  is 
simply  the  Green’s  function  (2.29)  or  (2.33).  The  first-order  term  in  the  perturbation 
series  is  found,  by  substituting  in  the  zeroth-order  solution,  to  be 

tfa  =  L[fa]  =  L\p0]  .  (2-45) 

The  second-order  term  is  similarly  found  to  be 

e2fa  =  L[e fa]  =  L2\p0]  .'  (2.46) 

Continuing  in  this  manner,  the  perturbation  series  for  the  total  field  is 

p  =  Po  +  L\po]  +  L2[po]  +  £3[po]  4 - •  (2.47) 

Redefining  the  scattered  field  as  a  series  of  higher-order  correction  terms, 

Ps  =  +  P(2)  4-  P(3)  4-  •  •  * 


? 


(2.48) 
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where  pU)  is  of  order  e,  is  of  order  e2,  and  so  on,  it  follows  that  each  higher-order 
correction  can  be  found  uniquely  in  terms  of  the  preceding  order  correction: 

p<n+1)  =  L\p^]  .  (2.49) 

Using  equation  (2.39),  the  perturbation  terms  for  the  volume  scattering  integral 
representation  are  found  to  be 

p(n+1)(r,*)  =  /  [P7*(r', t)g0(r, r>(n)(r', t)  + 

Jv  (2.50) 

•  7p(r',  t)Vg0(r,  r')  •  V'p(n)(r'>  t)]  dr' 

for  n  =  0, 1, 2,  •  •  • .  If  the  fluctuations  are  within  the  range  of  validity  of  the  pertur¬ 
bation  solution,  a  low-order  approximation  will  sufficiently  represent  the  scattered 
field. 

The  equivalent  perturbation  result  can  also  be  found  by  direct  iteration  of  the 
integral  equation  [67].  It  can  also  be  found  by  substituting  the  perturbation  series 
directly  into  the  wave  equation  for  the  LHS  and  satisfying  boundary  conditions  along 
the  interface  [42].  The  methods  are  equivalent  and  need  not  be  repeated. 

2.3.1  First-Order  Perturbation  Method 

For  a  source  and  receiver  positioned  in  the  UHS  and  in  the  far-field  of  the  interface 
(see  Figure  2.1),  the  first-order  scattered  field  can  be  found  in  terms  of  a  convenient 
Fourier  integral  over  the  scattering  volume.  Using  the  method  of  stationary  phase 
to  evaluate  the  Green’s  function  (2.29)  in  the  far  field  (see  Appendix  B),  the  half¬ 
space  Green’s  function  for  transmission  from  the  lower  to  the  upper  medium  is 
approximated  as 

flb(r,rO  =  — .  (2.51) 

p  47rr 

Similarly,  the  incident  field  at  r'  in  the  LHS  due  to  a  source  at  ro  is 

Po(  r')  =  <7o(r',  r0)  =  TMe^e^-r'  (2.52) 

47rro 
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The  position  vectors  r  and  r0  lie  in  the  water,  and  r'  lies  in  the  sediment.  The  wave 
vectors  k,  and  ks  are  oriented  in  the  incident  and  scattering  directions  in  the  LHS, 
respectively.  They  are  defined  as 


kj  =  kxix  +  kyiy  -  fyz  , 
k§  =  kxsx  d"  kysy  "f“  P$z  , 


(2.53) 


with  pi  and  ps  defined  analogously  to  (2.30).  The  horizontal  components  of  the 
incident  and  scattered  field  wavevectors  are  defined  as  kxi  =  kox0/r0,  kxs  =  k0x/r, 
kyi  =  k0y0/rQ,  and  kys  =  k0y/r.  The  coefficient  T  is  the  pressure  transmission 
coefficient  for  a  plane  wave  incident  upon  the  interface  from  the  water  at  the  incident 


and  scattering  angles,  0,  and  6S  respectively: 

2  pPo 


m )  = 

T(0S)  = 


PoPi  +  pPo 

2pPo 


(2.54) 


PoPs  +  pPo 

Using  (2.50),  (2.51)  and  (2.52),  the  first-order  bistatic  field  at  position  r  is  found 
to  be 


p">(r ,()  =  r(r)  (  dr'[e^(r',t)  +  7(,(r',()(k.- •  k,)] ,  (2.55) 

Jv 

where 

P('r\  _  Po  »fco(ro+r)  (2.56) 

U }  -  p  (4t r)2r0r  ’  V  ; 

Figure  2.3a  illustrates  the  first-order  result  graphically.  The  incident  field  is  refracted 

into  the  lower  medium  and  interacts  with  the  inhomogeneities.  The  scattered  field  is 

an  integral  over  all  the  monopole  (7*)  and  dipole  (jp)  sources  in  the  volume  shifted 

in  phase  by  the  factor  related  to  the  incident  and  scattered  directions.  The  integral 

results  in  a  Fourier  transform  of  the  medium  fluctuations  over  the  volume  of  the 

medium.  Since  the  volume  of  the  scattering  medium  is  the  half-space,  the  Fourier 

integral  must  be  windowed  in  the  depth  direction. 
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The  equivalent  result  for  a  two-dimensional  medium  is  found  using  the  two- 
dimensional  far-field  Green’s  functions 


9o(r,r') 


■  Po  T(0S)  ^ikorg-iks-r' e~in/4 

p  \ft£irkor 


(2.57) 


and  incident  field 


po(r')  =  -7-  T(6i}  ..,Pikor0fjkt.r'f>-^/4  ?  (2.58) 

y/8Trk0r0 

where  r  =  ( x,z )  (see  Appendix  A).  The  two-dimensional  integral  equation  for  the 
scattered  field  is  the  same  as  (2.55)  except  the  geometric  factor  is  redefined  as 

r(r)  =  ^T(^)T(^)etfc0(ro+r)  ^  (2.59) 

p  8Ttk0^F 

and  the  integration  is  over  the  two-dimensional  medium. 


2.3.2  Higher-Order  Perturbation  Fields 


The  scattered  fields  represented  by  higher-order  terms  in  the  perturbation  series  can 
be  interpreted  physically  by  examining  the  recursive  nature  of  (2.50).  The  second- 
order  perturbation  result,  which  contains  a  double  integral  operation,  is  found  by 
substituting  the  first-order  result  into  (2.50): 


P(2)(r)  =  k4  JJ 7K(r2)7*(ri)0o(r,  r2)#0(r2,  r1)g0{r1,r0)dr1dr2 

k2  JJ 7«(r2)7p(ri)2o(r,  r2)V1^o(r2,r1)  •  Vi^0(ri,r0)dr1dr2 

k2  JJ 7p(r2)7/c(ri)V2p0(r,  r2)  •  V2flf0(r2,  ri)^0(ri,r0)dridr2 

JJ 7p(r2)7p(ri)V2£f0(r,r2)  •  [V2Vi0O(r2,ri)]  •  Viflb(ri,r0)«fPi</r2 


+ 
+  . 
+ 


(2.60) 


The  operator  Vx  is  the  gradient  with  respect  to  iq;  V2  is  the  gradient  with  respect  to 
r2;  and  V2Vi  is  a  dyad  operator.  The  time  variable  t  is  suppressed  for  convenience. 
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Figure  2.3:  Small  perturbation  method  for  volume  scattering  from  a  random  half¬ 
space:  (a)  first-order  scattering;  (b)  second-order  scattering. 
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Consider  the  double  integrals  in  the  second-order  expression.  The  presence  of  the 
Green’s  function  <70(1*2, 1*1)  represents  the  interaction  between  the  elemental  volumes 
dr  1  and  dr2.  Therefore,  the  second-order  term  in  the  perturbation  series  represents 
the  field  due  to  second-order  scattering  in  the  medium,  that  is,  scattering  due  to  the 
interaction  of  two  elemental  volumes.  The  first  integral  is  double  scattering  due  to 
compressibility  fluctuations.  The  second  and  third  integrals  are  double  scattering 
between  compressibility  and  density  fluctuations.  The  fourth  integral  is  double  scat¬ 
tering  due  to  density  fluctuations.  The  dot  products  of  the  gradients  of  the  Green’s 
functions  can  be  thought  of  as  creating  a  dipole  directivity  of  the  scattered  field.  In 
the  far-field,  the  gradient  in  the  Green’s  function  creates  dipole  patterns  oriented  in 
the  direction  of  propagation.  Figure  2.3b  illustrates  the  second-order  interactions 
due  to  scattering  by  density  and  compressibility  fluctuations.  The  third-order  per¬ 
turbation  term  includes  third-order  scattering.  Similarly,  higher-order  terms  of  the 
perturbation  series  represent  higher-order  scattering  in  the  medium.  The  scattered 
field  (for  each  higher  order)  will  be  composed  of  all  the  possible  compressibility  and 
density  interactions. 

Scattering  diagrams  (due  to  Feynman)  are  a  useful  method  of  expressing  the 
higher-order  terms  [68,  80].  Let  <70(1*1,  r0)  be  represented  by  a  straight  line  segment: 

<7o(r2,  ri)  ~  - •  (2-61) 

r2  1*1 

The  vertices  of  the  line  segment  represent  the  integration  variables.  Several  new 
symbols  are  defined  to  represent  the  medium  and  derivatives  of  the  Green’s  func¬ 
tion.  A  filled  vertex  represents  a  compressibility  fluctuation,  and  a  hollow  vertex 
represents  a  density  fluctuation: 

£27«(r)  ~  •  , 

7 pOO  ~  °  • 

An  arrow  on  the  line  segment  represents  the  gradient  of  the  Green’s  function  with 
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respect  to  the  integration  variable  the  arrow  is  pointing  towards: 


Vi0O(r2,ri) 

V2firo(r2,ri) 


(2.63) 


The  second-order  perturbation  term  can  now  be  expressed  in  diagrams  as 
p(2)(r)  = - • - • - 1 - • — y-o-i - h 

r  1*2  ri  ro  r  *2  ri  ro 


(2.64) 


r  **2  ri  ro  r  **2  **i  ro 


The  diagrams  can  be  simplified  by  dropping  the  coordinates  of  the  inner  vertices. 
Noting  also  that  the  hollow  vertices  are  always  bracketed  by  arrows  representing  the 
gradient  operation,  the  arrows  will  be  dropped  and  the  derivatives  implied.  The 
third-order  perturbation  result  is  represented  as  follows: 


P(3,(  r)  = 


— +  - 
ro  •  r 

— +  - 
r0  r 

—  +  - 
r0  r 

—  +  - 
ro  r 


(2.65) 


2.4  Temporal  and  Spatial  Correlation  of  the  Scattered  Field 


As  the  sediment  fluctuates  randomly  in  space  and  time,  the  scattered  field  will 
fluctuate.  Thus,  the  scattered  field  is  a  random  process,  which  can  be  characterized 
by  the  cross-correlation  of  ps  at  two  different  locations  and  times: 

CPP  =  (p.(ri,<i)p*(r2,*2))  •  (2.66) 

Using  the  perturbation  result  (2.55)  and  considering  only  the  case  where  ri  =  r 2, 
the  correlation  of  the  first-order  bistatic  scattered  field  can  be  found.  Assuming 
the  randomness  in  the  medium  is  stationary  in  space  and  time  and  making  the 
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Zj 


Figure  2.4:  Half-space  windowing  in  the  difference  coordinates. 

change  of  integration  variables  to  the  difference  and  center  coordinates,  defined  as 
rd  =  r'x -r2,  and  r'c  =  (r'x  +r2)/2,  respectively,  and  the  difference  in  time  r  =  f i  -f2, 
the  correlation  function  is  found  to  be 

Cpp(t)  =  A|r(r)|2  [°  dz'c  f  dr'd  C (rd', r )e~ikd'rd' e~2az'c  .  (2.67) 

The  difference  between  the  incident  and  scattered  wave  vectors  is  the  Bragg  wave 
vector,  krf  =  Re[k,  —  kj.  The  area  of  the  scattering  region  is  A,  and  a  is  the 
imaginary  part  of  the  vertical  component  of  the  difference  between  the  incident  and 
scattered  wave  vectors: 

a  =  Im[(k,-  -  ks)  •  z]  =  -Im [/?,•  +  (da\  .  (2.68) 

The  integral  on  the  difference  coordinate  is  over  the  windowed  half-space  volume 
(see  Figure  2.4)  defined  as 

p  poo  poo  p—2zc 

/  drd=  dx'd  /  dyd  /  dzd  .  (2.69) 

J  hs  J —oo  J —oo  J  2zc 

The  combined  sediment  correlation  function  is 

C(rd,r)  =|£|4  CKK(rd,r)  +  2Re[P(kj  •  ks)]  CKp{rd,r )  + 

|kf  •  ks|2  Cpp(rd,T)  . 


(2.70) 
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The  functions  CKK,  CKp,  and  Cpp  are  the  auto-  and  cross-correlations  of  the  zero-mean 
functions  defined  in  Section  2.1.1. 

If  it  is  assumed  that  the  medium  compressibility  is  proportional  to  the  density 
using  (2.14),  then  the  correlation  function  is  simplified  to 

C( U,T)  =  \pk2  +  (k;  ■  ks)|2  Cpp{Vd,r)  .  (2.71) 

The  effect  of  the  half-space  geometry  on  the  correlation  function  (2.67)  is  to 
create  a  sliding  windowing  in  the  zd  direction,  as  defined  in  Figure  2.4.  The  half¬ 
space  windowing  function  is  a  depth-dependent  rectangular  window, 


w(zd,  zc)  = 


having  the  Fourier  transform 


W(kz,zc) 


1  if  \z<i\  <  2zc, 

0  otherwise  , 


sin  (2  kzzc) 


(2.72) 


(2.73) 


Applying  the  Wiener-Khintchin  theorem  and  the  Fourier  windowing  theorem, 
the  correlation  function  (2.67)  can  be  expressed  in  terms  of  the  time-domain  spatial- 
spectrum  of  density  fluctuations: 


Cpp(r)  =  8t t3A  |r(r)|2  \fik2  +  (kt-  •  ks)|2  Spp( kd,r)  . 


(2.74) 


The  density  spectrum,  Spp,  is  not  the  true  spectrum  but  rather  the  convolution 
of  the  true  density  spectrum  and  the  Fourier  transform  of  the  half-space  window, 
defined  as 


/0  poo 

dzc  /  dk'dz  Spp(kdx,  kdy,  kdz,  r)W{k’dz  -  kdz,  zc)e~2aZc  .  (2.75) 

•oo  J —oo 

The  convolution  integration  of  (2.75)  is  only  over  the  ^-component  of  the  Bragg 
wavenumber.  Using  (2.73),  the  integration  over  zc  in  (2.75)  can  be  performed  to 
find  the  modified  spectrum: 


Spp  (krf,  t 


Sppi kdx  1  kdy  1  kdz,  t)  j 

2*[«2  +  (k'dz-kdzf]  * 


(2.76) 
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The  modified  spectrum  can  also  be  shown  to  be  the  Fourier  transform  of  the  true 
density  correlation  function  windowed  in  depth  by  the  medium  attenuation: 

s„,( ki,T)  =  •  (2-77) 

Therefore,  the  half-space  effectively  modifies  the  correlation  structure  of  the  medium 
in  the  depth  direction.  If  attenuation  in  the  medium  is  low,  the  half-space  effects 
will  be  minimized.  The  incident  field  will  penetrate  a  significant  depth  into  the 
medium,  relative  to  the  correlation  length  of  the  fluctuations,  and  scattering  from 
the  half-space  medium  can  be  modeled  using  the  true  spectrum.  If  the  attenuation 
is  high,  the  field  will  be  limited  to  a  layer  near  the  interface.  If  the  penetration 
depth  is  smaller  or  comparable  in  dimension  to  the  correlation  scale  of  the  medium, 
the  half-space  effects  will  be  significant. 


2.5  Scattering  Cross-Section 

The  scattering  cross-section  is  defined  in  terms  of  the  incoherent  intensity  of  the 
scattered  far-field.  If  the  penetration  distance  into  the  sediment  is  small  compared 
to  the  distance  traveled  in  water,  it  is  typical  to  describe  volume  scattering  in  terms 
of  an  equivalent  surface  scattering  strength.  For  a  pressure  wave  incident  on  the 
surface  of  area  A ,  the  bistatic  surface  scattering  cross-section  (per  unit  area  per  unit 
solid  angle)  is  defined  as 


Is{r)r2 

<Ts  =  ~Ta~' 

where  r  is  the  distance  to  the  receiver  in  the  far-field  (Figure  2.1),  and 


(2.78) 


(2.79) 


is  the  incident  intensity  immediately  above  the  interface.  For  convenience,  crs  will 
simply  be  referred  to  as  the  cross-section.  The  scattering  strength  is  defined  as 
101og10(crs)  in  units  of  dB. 
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The  incoherent  scattered  intensity  is  defined  in  terms  of  the  zero-lag  covariance 
of  the  scattered  field  as 

/. M  =  ^  [(P.MpJM)  -  l(p.M)|2]  •  (2.80) 

The  scattered  field  covariance  can  be  found  using  the  perturbation  expansion  for 
the  scattered  field  (2.48).  By  consistently  grouping  terms  in  the  resulting  covariance 
expansion  by  powers  of  the  small  parameter,  the  scattering  cross-section  is  found  as 
the  perturbation  series 

ers  =  a*1)  +  erj2)  +  cr<3)  +  cr^  +  •  •  •  .  (2.81) 

Each  higher-order  term  in  the  cross-section  series  corresponds  to  the  complete  group¬ 
ing  of  like  orders  from  the  scattered  field  covariance.  The  first-order  cross-section  is 
zero.  The  second-order  cross-section  is  found  to  be 

^2)  =  3£p[<|P(,,|2)-»S2)]  .  (2-82) 

which  is  second-order  in  the  scattered  field,  and 

</i2)  =  I(P(1)>|2  •  (2-83) 

The  third-order  cross-section  is 

[<P(  V(2)  +  Pl  V(1)>  -  Vf  >]  .  (2-84) 

where 

<ti3)  =  (p(i)Xp'(2)>  +  (P(2)}<P‘W>  ■  (2-85) 

The  fourth-order  cross-section  is 

<T*4>  =  [dP<2)l2  +p'V(3)  +P(V(1))  -  fi4’]  ,  (2.86) 
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where 


if  =  l(p(2))l2  +  (p(1)Kp‘(3)>  +  (p(3)Xp*(1,>  • 


(2.87) 


Note  that  ( p f1))  =  0,  and  for  a  Gaussian  process  all  odd-order  cross-sections  are  also 
zero. 

Using  (2.55)  for  the  first-order  scattered  field,  and  using  (2.52)  to  find  the  inci¬ 
dent  intensity,  the  second-order  surface  scattering  cross-section  as  a  function  of  the 
incident  and  scattering  angles  is  found  to  be 


*!2)(m.) = m  \mw  im)i2  + (k,  •  k.)i2sw(kd) . 

L  p 

Defining  a  depth-dependent  volume  cross-section  as 


0-5  =  /  <Fv{Zc) 

J  —OO 


(2.88) 


(2.89) 


and  using  (2.75),  the  second-order  volume  cross-section  is  found  as 

of  (*)  =  5I/1P  +  (k.  •  k,)|2  r  S„(K)W(K  -  *.,*«)  K  ■  (2-90) 

^  J- OO 

This  expression  for  cross-section  is  similar  to  those  derived  by  [32,  34,  26],  except 
half-space  effects  are  included. 

The  total  scattering  cross-section  of  the  volume  is  defined  as 


tv  -  /  &v 

J  4tt 


(2.91) 


where  dui  is  the  differential  solid  angle.  The  total  surface  scattering  cross-section  is 


defined  as 


=  f  <r» 

J  2?r 


(2.92) 


where  du>  is  the  differential  solid  angle  over  the  UHS  hemisphere. 

In  two  dimensions,  the  second-order,  bistatic,  surface-scattering  cross-section  is 
defined  as 


vs(0i,0s)  = 


IiL  ’ 


(2.93) 
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where  L  is  the  length  of  the  horizontal  scattering  region.  The  two-dimensional  first- 
order  scattering  cross-section,  using  (2.57)  and  (2.58),  is  found  to  be 

<^(0.,  ».)  =1-2  T  lr(0i)|2  |r(0.)|2  W  +  (ki  ■  k,)|2  $„( kd,  r)  ,  (2.94) 

Z  p  rC  o 

where  Spp  is  now  the  two-dimensional  half-space  spectrum. 


2.6  Scattering  of  a  Narrow  Band  Pulse 


In  many  practical  applications  of  scattering  from  a  random  medium,  it  is  necessary 
to  use  a  finite  duration  pulse  as  the  incident  field  and  consider  the  time-domain  scat¬ 
tered  field.  The  scattered  field  is  typically  time-gated  to  correspond  to  scattering 
from  an  isolated  and  finite  region  of  the  scattering  medium.  This  is  especially  true 
for  scattering  from  the  seafloor  where  multipath  propagation  can  create  complicated 
interactions  of  different  scattering  mechanisms.  It  is  necessary,  therefore,  to  relate 
the  time-independent  scattering  of  a  continuous  wave  field  to  the  time-dependent 
scattering  of  a  pulsed  field.  In  this  section,  the  classical  narrow-band  approxima¬ 
tion  for  a  time-varying  medium  [30]  is  applied  to  the  half-space  volume  scattering 
problem. 

Consider  the  scattered  field  ps(t)  at  a  point  in  space  due  to  a  pulsed  input  field 
Pi(t),  suppressing  the  spatial  dependence  of  the  field  for  a  moment.  The  fields  are 
defined  as  the  real  component  of  the  complex  envelope  functions  modulating  the 
carrier  at  frequency  ujq: 

ps(t)  =  Re  [us(t)e~'“ot]  ,  ^295) 

Pi(t)  =  Re  [ui(t)e~tWot]  . 

The  envelope  functions  can  be  expressed  in  terms  of  a  time-varying  amplitude  and 
phase: 


u3(t)  =  As(t)e^  , 
Ui(t)  =  . 


(2.96) 
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For  a  narrow  band  input  signal,  AJt)  and  <f>i(t)  are  slowly  varying  functions  of  time 
relative  to  the  carrier  frequency. 

Since  the  scattering  medium  is  linear,  the  scattered  output  field  can  be  related 
to  the  input  field  by  the  integral 


Ps(t)  =  [  Pi{t')h(t,t')dt'  , 
J  —CO 


(2.97) 


where  h(t,  t')  is  the  response  of  the  medium  at  time  t  due  to  an  impulse  input  at  t' . 
Equivalently,  the  output  envelope  function  can  be  expressed  in  terms  of  the  Fourier 
transform  of  the  input  envelope  function  and  the  system  frequency  response: 


where 


/OO 

Ui(uL>)H{u o  +  u>,  t)e~tutdu> 

•OO 

1  f°° 

Ui(u)  =  —  j  Ui(t')e'wt' dt' , 


(2.98) 


(2.99) 


1  f°° 

H(u,t)  =  —  /  h(t,t  —  t')eluJt' dt' 

Jo 


(2.100) 


The  function  H (u>,  t)  is  the  time-varying  transfer  function  of  the  medium  (see  [30] 
Chapter  6).  By  definition,  it  is  the  scattered  output  field  when  the  input  field  is  a 
time-harmonic  function.  Therefore,  H(u,t)  is  simply  the  CW  solution  (found  using 
perturbation  or  exact  methods,  for  example).  Using  the  first-order  perturbation 
result  from  Section  2.3.1,  the  response  of  the  volume  scattering  medium  is 

H{u>,  t )  =  r(r)  /  [P7«(r#,  t)  +  7,(r',  f )(k  •  k.)]  e^-^'dr'  .  (2.101) 

Jv 


If  the  medium  is  time-invariant,  H(oj,  t)  is  time-invariant.  For  the  cases  of  interest  in 
this  thesis,  the  medium  is  varying  in  time  much  slower  than  the  acoustic  frequency, 
and  H(u),t)  can  be  considered  a  constant  in  time  over  the  duration  of  the  pulse 
scattering. 
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As  the  input  pulse  propagates  through  the  scattering  medium  and  interacts  with 
spatial  inhomogeneities,  the  scattered  output  field  becomes  a  random  function  of 
time.  Temporal  statistics  of  the  output  field  reflect  the  temporal  and  spatial  statis¬ 
tics  of  the  medium  and  the  characteristics  of  the  input  pulse.  The  pulse  character¬ 
istics  include  the  receiver  and  transmitter  beam  patterns,  the  pulse  duration  and 
bandwidth.  For  backscattering  at  any  instant  in  time,  the  region  enveloped  by  half 
(round  trip)  the  spatial  extent  of  the  pulse  becomes  the  effective  scattering  volume 
that  contributes  to  the  scattered  field  at  that  instant  in  time.  This  region  is  called 
the  “ensonified  volume.” 

A  fundamental  statistical  descriptor  that  shows  the  influence  of  the  ensonified 
volume  is  the  correlation  of  the  complex  envelope  of  the  scattered  field: 

Cuu(ti,t2)  =  (us(fi)u*(f2))  •  (2.102) 

In  this  case,  the  scattered  field  is  the  response  of  the  medium  to  a  single  pulse, 
and  the  times  t\  and  £3  both  occur  during  the  scattering  from  the  same  input  pulse. 
Between  the  correlation  times  t\  and  £2?  the  medium  is  assumed  to  be  time-invariant. 
Using  equation  (2.98)  with  the  time-dependence  of  the  medium  removed,  the  output 
correlation  is 

/oo  poo 

dui  /  du2  Ui(u>i) U* (u>2 ) CHh{ui , a>2 ) e-^1  +i“2 *2  ,  (2.103) 

■co  J  —oo 

where  Chh  is  the  two-frequency  correlation  function  of  the  medium,  defined  as 

Chh{w  1,^2)  =  (i/(w 0  +^i)^*(tt;o  +W2))  •  (2.104) 

For  volume  scattering,  using  (2.101),  and  applying  the  narrow  bandwidth  approxi¬ 
mation  for  the  input  pulse,  such  that 


(u>0  +  U>)2  «  OJq  , 


(2.105) 
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the  two-frequency  correlation  function  can  be  approximated  as 

dvd  x 

(7(r<i)  e~2aZc  e~iRe^'Vd  e-t*Cdl'(rc+rd/2)  e -*kd2-(rc+rrf/2)  ^ 


Cj hh(ui ,  u;2 )  =47r2  A  | T(r)  1 2e-1(u'1  -^)(r0-r)/c0  f°  ^  f 

J  —  OO  t/  Zl5 


(2.106) 


where  rc  is  the  location  of  the  center  of  a  patch  of  area  A.  Note  that  A  is  a  small 
area  (defined  by  performing  the  xc  and  yc  integrations  in  (2.106)),  not  the  entire 
ensonified  area  of  the  seafloor  as  the  pulse  travels  in  time  [47].  As  in  section  2.4, 
the  function  C(rd )  is  the  correlation  function  of  the  stationary  random  fluctuations 
in  the  medium.  The  wave  vector  k<£  is  defined  in  terms  of  the  carrier  frequency  u>o- 
The  wave  vectors  k^i  and  k <*2  are  defined  in  terms  of  and  u>2,  respectively. 

The  correlation  of  the  output  signal  is  then  found  using  (2.103)  and  (2.106)  and 
performing  the  integrations  with  respect  to  wi  and 

Cuuitnh)  =47rM|r(r)|2  /  dzc  [  drdx 

J— oo  J  hs  (2.107) 

Ui(t  1  -  -  t")  C( rd)  e-^ce-i^\kd].rd  _ 


The  time  delays  t'  and  t"  are  functions  of  the  distance  traveled  through  the  UHS 
and  into  the  LHS  medium, 


t'  =  ( i  -  s)  •  (rc  +  rrf/2)/e  -  t0 
t"  =  (i  -  s)  •  (rc  -  rd/2)/c  -  t0  , 


(2.108) 


where 


t0  =  (r0  -  r)/co  . 


(2.109) 


The  delay  times  are  complex  due  to  attenuation  in  the  LHS.  The  unit  vectors  i  and 
s  define  the  incident  and  scattering  directions  in  the  medium,  such  that 


(2.110) 
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Figure  2.5:  Scattering  of  a  pulse  of  duration  T0  from  a  random  half-space. 


The  translated  and  distorted  pulse  replica  Ui(t  —  t')  in  (2.107)  is  defined  using  an 
inverse  Fourier  transform  with  a  complex  time  delay: 

/OO 

.  (2.111) 

■OO 

In  general,  equation  (2.107)  can  be  considered  a  convolution  of  the  ensonified 
volume  and  the  medium  correlation  function  over  the  volume  of  the  scattering  region. 
The  pulse  volume  can  be  complicated  due  to  refraction  across  the  interface  and 
attenuation  in  the  lower  medium.  However,  if  the  pulse  envelope  is  a  simple  square 
wave,  the  ensonified  volume  can  be  defined  approximately  in  terms  of  its  pulse 
duration  T0  (see  Figure  2.5).  If  the  ensonified  volume  in  the  sediment  is  large 
compared  to  the  correlation  scale  of  the  medium, 

(2.112) 

the  convolution  can  be  ignored.  The  correlation  function  for  the  time-domain  scat¬ 
tered  field  can  then  be  approximated  using  the  single  frequency  correlation  result  of 
(2.67). 
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2.7  Volume  Scattering  from  a  Periodic  Random  Medium 

It  is  sometimes  convenient  to  assume  a  plane  incident  wave,  as  opposed  to  a  spher¬ 
ical  wave  originating  at  a  point  source  such  as  (2.52).  This  will  become  especially 
apparent  in  Chapter  4  when  numerical  methods  for  volume  scattering  will  be  dis¬ 
cussed.  However,  by  assuming  an  incident  plane  wave,  the  horizontal  extent  of  the 
scattering  volume  must  be  extended  to  infinity.  A  simple  method  of  handling  the 
infinite  horizontal  dimension  is  to  assume  periodicity  of  the  random  medium  in  the 
horizontal.  This  is  fundamentally  a  different  problem  than  volume  scattering  from 
the  seafloor,  since  the  seafloor  is  not  periodic.  However,  by  assuming  the  period  is 
much  larger  than  the  correlation  length  of  the  medium,  a  single  period  will  contain 
a  sufficient  statistical  description  of  the  medium.  Then,  the  scattered  field  from  the 
periodic  medium  can  be  assumed  to  be  representative  of  a  nonperiodic  medium  with 
the  same  statistics.  This  is  only  true  at  particular  scattering  angles  (corresponding 
to  the  Floquet  modes),  as  will  be  discussed. 

Therefore,  consider  a  plane  wave  incident  on  an  infinite  two-dimensional  half¬ 
space.  If  the  lower  half-space  is  periodic  in  the  ^-direction  with  period  L,  such  that 
the  random  fluctuations  in  the  medium  satisfy  the  conditions 

'yK(x,  z )  =  +  f,  2)  , 

7 p(x,z)  =  7 P(x  +  L,z)  , 

and  lc  <  L,  then  using  Floquet’s  theorem  [29,  80]  the  scattered  field  can  be  written 
as 

(2.114) 


p(x,  z)  =  p(x,  z)etkixX 


(2.113) 


where  p(x ,  z)  =  p(x  +  L,  z )  and  k{x  =  ko  cos  The  integral  equation  for  scattering 
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from  the  periodic  volume,  similar  to  (2.38),  is 


pL/2+mL  rO 

{x,z)=  ^  /  dx  dz'[P'jK(x,,z')g0(x,z-,x',z,)p(x',z') 

L/2+mjD  J — oo 


-L/2-\-mL  J — oo 


(2.115) 


+  7 'p(x',  z')Vg0(x,  z;  x',  z')  •  Vp(x',  z')]  . 


With  the  change  of  variables  x"  —  x'  —  mL  and  z"  =  z',  and  applying  the  periodic 
boundary  conditions,  the  integral  can  be  rewritten  as 


/L/2  [0 

dx"  /  dz"  [k2iK(x",  z")g0(x,  z;x",  z")p(x",  z") 

-L/2  J — oo 

+  i,(x\  z")V"go(x,  z I  x",  z")  ■  V'p(x",  *")]  , 
where  go  is  the  periodic  Green’s  function  defined  as 

OO 

g0(x,  z;  x",  z")  =  ^  g0(x,z-,x" +  mL,z")etkixmL  . 


(2.116) 


(2.117) 


Using  the  definition  of  the  half-space  Green’s  function  in  two  dimensions,  similar  to 
(2.33),  and  applying  the  Possion  sum  formula, 


e^-k^mL  =  -f  5^  *(**  “  ki*  ~  2*m/L)  , 


(2.118) 


the  periodic  half-space  Green’s  function  for  the  LHS  is  found  to  be 
g0(x,z;x",z")  =  -^-  N+-"l  +  ,  (2.119) 


where 


^ xm  —  ^ ix  ”1”  C2i'K7TI  j  Ij 


(2.120) 


Pm  —  \/k2  k%Ti 


(2.121) 


correspond  to  the  Floquet  modes  of  the  periodic  medium. 
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The  periodic  Green’s  function  for  transmission  from  the  lower  to  the  upper 
medium,  similar  to  (2.29),  is  found  to  be 


OO  qnf 


g0(x,  z]  x",  z") 


m  ^  A 

v  '  rar-rn 


T_  (hem)  eiPomz-i0mz" eikxm(x-x")  ^  (2.122) 


where  p0m  =  \Ao  "  klm- 


2. 7. 1  Scattering  Cross-Section  of  a  Periodic  Random  Medium 

In  the  strict  sense  of  the  definition,  the  scattering  strength  does  not  exist  for  an 
infinite  periodic  medium  because  the  far  field  of  an  infinite  scattering  region  does 
not  exist.  However,  in  a  practical  sense,  defining  a  scattering  strength  is  still  useful, 
especially  for  comparisons  with  non-periodic  scattering.  Therefore,  assume  the  pe¬ 
riodic  medium  is  of  some  finite  and  large  extent.  The  total  field  scattered  from  the 
finite  periodic  half-space  can  then  be  approximated  as  the  sum  of  all  the  propagating 
Floquet  modes  of  an  infinite  medium, 

n2 

p,(x,z)  =  £  ,  (2.123) 

m=Ni 

where  Ni  and  N2  define  the  range  of  modes  for  which  (30m  is  real.  The  scattered 
field  for  each  mode  pam  is  defined  using  (2.116)  and  the  Green’s  function  (2.122)  as 

rL!  2  fo 

Psm=  dx'  dz'  [fc27«(®',  z')gm {x',  z')p{x\  z') 

J-L/2  J-oo  (2.124) 

+  lp(x',z’)V'gm(x',z')-Vp(x',z')}  , 

with 


/V  _  j  Hkm)  p-ipmz'-ikxmx' 

l(  ’  j  2  L  An 


(2.125) 


The  differential  scattering  cross-section  as  a  function  of  scattering  angle  is  found 
by  relating  the  total  scattered  field  (2.123)  to  the  total  scattering  cross-section. 
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The  total  equivalent  surface  scattering  cross-section  for  the  region  is  defined  as  the 
integral  of  the  differential  cross-section  over  the  scattering  directions, 


=  J  as(Os)dOs 

=  fko  (Ts{kx) 

J_kok0sm6s  *  ’ 


(2.126) 


with  kx  =  k0  cos  6S.  For  the  periodic  medium,  the  integral  is  replaced  by  a  summation 
over  the  propagating  modes, 


vA  ^sjkxm)  2?r 
ko  sin  0sm  L 

m=Ni 


(2,127) 


with  dkx  =  27 t/L,  and  with  0sm  being  the  scattering  angles  in  the  far-field  associated 
with  the  propagating  modes. 

Equivalently,  the  total  scattering  cross-section  can  be  defined  in  terms  of  the 
incident  and  scattered  power  flux  density  [29]  as 


fsn  • da 


(2.128) 


where  the  area  is  an  arbitrary  surface  enclosing  the  UHS,  and  da  is  the  differential 
area  vector  normal  to  the  surface.  The  area  integral  will  be  rewritten  as  an  integral 
of  the  normal  component  of  the  power  flux  over  a  line  just  above  the  interface, 


CTts  — 


IZo  h-zdx 


(2.129) 


where 


Using  the  plane  wave  intensity  of  the  incident  field, 


(2.130) 


2poCo  ’ 


(2.131) 
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and  exploiting  the  orthogonality  of  the  Floquet  modes  in  (2.123)  to  perform  the 
integration  over  x ,  the  total  scattering  cross-section  is  found  in  terms  of  a  summation 
over  propagating  modes: 

JV2 

c Tts=  (|psm|2)sin0sm  .  (2.132) 

m=Ni 

Equating  expressions  (2.132)  and  (2.127),  it  is  evident  that  the  discrete  differen¬ 
tial  cross-section  for  a  periodic  medium  can  be  found  as  a  function  of  the  discrete 
scattering  directions  associated  with  each  Floquet  mode: 

&s  =  7J  A:o  sin  $sm(|p«m|  )  • 

27 r 


(2.133) 
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Chapter  3 

NUMERICAL  METHODS  FOR  VOLUME  SCATTERING 

The  basic  solution  strategy  of  most  numerical  methods  is  to  reduce  the  equa¬ 
tions  governing  a  particular  problem  into  a  matrix  equation,  and  then  solve  for  the 
unknowns  using  efficient  matrix  inversion  techniques.  Such  methods  are  widely  ap¬ 
plied  to  wave  propagation  and  scattering  problems:  finite  difference,  finite  element, 
finite  volume,  etc.  Each  method  begins  by  formulating  the  wave  solution  as  either  a 
partial  differential  equation  with  specified  boundary  and  initial  conditions  or  by  for¬ 
mulating  an  integral  equation  solution.  In  this  chapter,  the  discussion  of  numerical 
solutions  to  the  volume  scattering  problem  will  be  limited  to  the  single  frequency 
integral  equation  solution  using  the  method  of  moments.  Overviews  of  this  and  other 
numerical  methods  for  wave  scattering  can  be  found  in  texts  on  the  subject  [24,  57] 
and  in  review  articles  with  extensive  bibliographies  [18,  81],  for  example.  For  an 
application  of  a  finite-difference  time-domain  solution  to  surface  scattering,  see  [25], 
and  for  finite- volume  methods  for  volume  scattering  see  [44]. 

The  selection  of  a  particular  numerical  method  should  follow  from  the  geometry 
and  boundary  conditions  of  the  problem  at  hand  and  the  type  of  solution  desired 
(time  or  frequency  domain,  far-field,  etc.).  Without  presenting  a  review  of  merits  and 
weaknesses  of  other  common  numerical  methods,  the  method  of  moments  (MoM)  is 
chosen  for  several  reasons.  For  one,  sediment  volume  scattering  is  generally  confined 
to  a  layer  near  the  sediment-water  interface  (due  to  loss  in  the  medium).  There¬ 
fore,  the  size  of  the  scattering  volume  and  the  resulting  number  of  unknowns  in  the 
numerical  solution  is  limited.  This  allows  a  MoM  solution  of  the  volume  scattering 
integral  equation  to  be  used.  In  general,  a  numerical  solution  using  an  integral  equa- 
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tion  formulation  would  not  be  well-suited  for  volume  scattering  problems  because  of 
the  size  of  the  solution  domain. 

Since  it  can  be  used,  the  MoM  solution  benefits  from  the  integral  equation  for¬ 
mulation  of  the  problem.  The  integral  operator  contains  the  Green’s  function  and 
therefore  the  complete  boundary  conditions  of  the  problem.  There  is  no  need  to 
approximate  the  discontinuity  and  boundary  conditions  at  the  half-space  interface 
by  discretization.  The  half-space  Green’s  function  is  calculated  separately  to  any 
desired  accuracy  using  standard  numerical  wavenumber  integration  techniques.  This 
eliminates  the  problems  of  grid  dispersion  suffered  by  FD  and  FEM  methods  [57]. 
The  use  of  the  Green’s  function  in  the  CW  integral  equation  also  properly  accounts 
for  absorption  in  the  medium.  Absorption  is  more  difficult  to  model  in  the  time 
domain.  Finally,  the  integral  equation  MoM  solution  only  requires  the  discretiza¬ 
tion  of  the  scattering  volume  rather  than  the  entire  volume  enclosing  the  scattering 
medium. 

Perhaps  the  most  important  benefit,  for  purposes  of  this  thesis,  is  that  the  MoM 
solution  of  scattering  using  the  IE  formulation  is  easily  compared  with  the  small 
perturbation  solution.  The  small  perturbation  method  for  volume  scattering  is  also 
formulated  as  an  integral  equation  (see  Chapter  2),  which  can  then  be  handled 
numerically  using  the  MoM  strategy.  The  resulting  matrix  operator  is  identical  to 
the  matrix  operator  of  the  exact  solution.  The  perturbation  solution  is  found  to 
various  orders  by  iteration,  while  the  exact  solution  is  found  by  inversion.  Thus, 

r 

the  MoM  solution  is  ideal  for  studying  the  accuracy  and  validity  of  the  perturbation 
method. 
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3.1  Method  of  Moments 

To  formulate  a  numerical  solution  for  volume  scattering,  the  integral  equation  (2.37) 
is  rewritten  in  operator  notation, 

p(r)  =  Po(r)  +  Lp{ r)  ,  (3.1) 

where  L  is  the  linear  integro-differential  operator  of  equation  (2.38).  The  operator 
L  has  a  weakly  singular  kernel  due  to  the  Green’s  function.  The  unknown  is  the 
field  variable  p(r).  Assume  that  a  solution  exists  and  is  given  by 

p(r)  =  (I-L)~1p0(r).  (3.2) 

The  basic  solution  strategy  is  to  reduce  the  linear  system  into  a  matrix  equation, 
and  then  solve  for  the  unknowns  using  efficient  matrix  inversion  techniques.  When 
the  field  position  r  on  both  sides  of  the  equality  is  inside  the  scattering  volume, 
(3.1)  represents  an  inhomogeneous  Fredholm  integral  equation  of  the  second  kind. 
In  general,  solution  techniques  are  stable  and  the  resulting  matrix  equations  well- 
conditioned.  In  addition,  the  kernel  of  the  integral  operator  containing  the  Green’s 
function  (2.33)  is  a  smooth  function  that  decays  with  distance  due  to  attenuation 
in  the  medium.  These  properties  provide  a  solution  matrix  which  in  general  will  be 
diagonally  dominant  and  potentially  sparse.  Efficient  numerical  methods  for  matrix 
operations  and  storage  can  then  be  exploited.  After  solving  for  the  unknown  field 
inside  the  volume,  the  field  outside  is  found  by  evaluating  the  integral  equation 
directly  using  quadrature  methods,  for  example. 

The  numerical  solution  begins  by  expanding  the  unknown  field  .in  a  series  of 
suitably  chosen  basis  functions  fn: 

OO 

Kr)  =  ■ 

n—  1 


(3.3) 
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The  function  set  {/„}  forms  a  true  basis  if  the  inner  product  space  forms  an  ( 

onal  set, 

>rthog- 

(fmJn)  =  n^m, 

(3.4) 

where  ()  is  the  inner  product  integral,  and  if  the  basis  set  satisfies  completeness. 
Completeness  in  this  context  means  that  any  well-behaved  function  can  be  repre¬ 
sented  by  the  series  to  any  degree  of  accuracy  [24],  or 

CO 

n=l 

(3.5) 

Approximate  representations  of  the  field  are  found  by  choosing  a  finite-dimensional 

subspace  of  the  basis  functions, 

N 

Mr)  =  ^Pnfnir)  , 

71=1 

(3.6) 

and  solving  for  the  unknown  coefficients  ( pn )  by  minimizing  the  residual 

N 

R(r)  =  p0{r)  -  ]Tpn[/n(r)  “  LMr)]  ■ 

71=1 

(3.7) 

Using  the  method  of  weighted  residuals  and  choosing  a  set  of  weighting  or  testing 

functions  wm ,  the  unknown  coefficients  can  be  found  by  forcing  the  residuals  to  be 
orthogonal  to  the  testing  functions  with  the  inner  product  integral  [57,  24], 

[  wm(r)R(r)dr  =  0,  m  =  1, 2,  •  •  • M  . 

Jv 

(3.8) 

The  testing  functions  do  not  necessarily  take  the  same  form  as  the  basis  functions. 

The  problem  now  reduces  to  the  solution  of  a  system  of  linear  equations,  and 

(3.8)  can  be  rewritten  as  the  matrix  equation 

[A  -  Z]p  =  p0  • 

(3.9) 
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The  matrices  Z  =  [ Zmn ]  and  A  =  [Amn]  are  of  size  M  x  N  with  elements 

Zmn  =  /  wm{r)Lfn{r)dr  , 

Jv 

A mn  =  J  wm(r)fn(r)dr  . 

The  vector  p0  =  [po,m]  has  elements 


(3.10) 


po,m  =  /  wm(r)po(r)dr  .  (3.11) 

Jv 

The  solution  vector  p  =  \pn]  is  comprised  of  the  unknown  coefficients  of  the  basis 
function  expansion.  The  matrix  Z  will  be  referred  to  as  the  scattering  matrix.  The 
major  complexity  in  finding  the  solution  now  lies  in  the  calculation  of  these  so  called 
moment  integrals  (hence  the  name  of  the  solution  method)  and  the  inversion  of  the 
resulting  matrix  A  —  Z.  The  complexity  involved  is  a  direct  consequence  of  the  choice 
of  basis  and  testing  functions. 


3.1.1  Perturbation  Solution 

The  small  perturbation  solution  of  Section  2.3  can  also  be  written  in  terms  of  a  ma¬ 
trix  equation.  Rewrite  the  perturbation  integral  equation  (2.50)  in  operator  notation 
as 

p(‘+1)(r)  =  Lp^\r)  ,  (3.12) 

where  L  is  the  same  operator  as  above.  Each  lower-order  term  in  the  perturbation 
series  is  then  expanded  as  a  series  of  basis  functions: 

p(i)«  S  £>£Vn(r)  .  (3.13) 

n= 1 

The  solution  for  each  higher-order  term  in  the  perturbation  series  involves  the  mul¬ 
tiplication  of  the  scattering  matrix  and  the  unknown  coefficients  of  the  preceding 
order: 


pb+i)  _  i .  pW  # 


(3.14) 
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The  solution  vector  p(l+1)  =  [pm+1']  has  elements 

Pm+1)  =  /  wm{r)p{i+1)(r)dr  .  (3.15) 

Jv 

The  scattered  field  for  the  perturbation  solution  can  then  be  found  by  iteration, 
where 

Ps  =  P(1)  +  p(2)  +  •  •  •  •  (3.16) 

Using  (3.14)  with  M  =  N,  the  total  scattered  field  field  is  found  in  terms  of  the 
scattering  matrix: 

OO 

p  =  ^Zn-p0.  (3.17) 

71=0 

For  a  weakly  fluctuating  medium,  the  series  will  converge  rapidly  and  can  be  trun¬ 
cated  to  approximate  the  complete  solution  found  from  (3.9). 

3.2  Discretization  and  Selection  of  Basis/Testing  Functions 

The  discretization  of  the  volume  scattering  integral  operator  into  a  finite-dimensional 
subspace  of  basis  functions  requires  approximations  that  compromise  the  solution 
accuracy  for  the  purpose  of  computational  efficiency.  The  choice  of  basis  functions  is 
the  principal  issue  in  implementing  the  MoM  solution,  as  discussed  in  [57]  and  [24]. 
The  selection  influences  the  accuracy  of  the  approximate  solution,  the  complexity  of 
the  computation  of  the  matrix  elements,  and  the  size  of  the  scattering  matrix.  While 
the  solution  with  the  best  accuracy  is  always  desired,  computational  limitations 
may  motivate  a  less  optimal  formulation.  In  some  cases,  computational  resources 
(memory)  put  an  upper  limit  on  the  matrix  size.  In  other  cases,  the  time  required 
to  perform  computations  that  create  the  scattering  matrix  limits  the  problem. 

The  choice  of  basis  and  testing  functions  is,  in  general,  a  function  of  the  geom¬ 
etry  and  boundary  conditions  of  the  problem,  the  expected  form  of  the  solution, 
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and  resources  at  hand.  Formally,  the  basis  functions  should  satisfy  conditions  of 
orthogonality  and  completeness,  as  discussed  previously.  However,  approximate  so¬ 
lutions  can  still  be  found  regardless  of  whether  wn  and  fn  form  true  basis  sets.  In 
practice,  the  functions  are  often  not  orthogonal.  This  creates  difficulties  in  deriv¬ 
ing  analytic  statements  about  the  numerical  convergence  of  the  approximation  [57]. 
Convergence  is  usually  shown  numerically  by  comparing  the  approximate  solution 
with  exact  solutions  as  N  ->  oo  . 

Many  types  of  basis  functions  are  applied  in  the  current  literature:  Dirac  delta 
functions,  rectangular  pulse  functions,  triangle  or  roof-top  functions,  quadratic  splines, 
and  Lagrangian  polynomials,  for  example.  Some  of  these  form  orthogonal  bases  and 
others  do  not.  In  general,  lower  interpolation  error  is  achieved  with  higher-order 
basis  functions  (e.g.,  high  polynomial  order).  However,  in  practice  a  trade-off  must 
be  made  between  accuracy,  computational  efficiency,  and  general  applicability  of  the 
solution  method  to  a  range  of  problem  geometries.  For  volume  scattering,  this  trade¬ 
off  is  especially  noticeable.  The  number  of  unknowns  required  to  provide  sufficient 
accuracy  can  become  very  large  (compared  to  similar  numerical  methods  used  to 
solve  surface  scattering  problems) .  Alternatively,  reducing  the  number  of  unknowns 
by  selecting  higher-order  or  full-domain  basis  functions  can  create  numerical  com¬ 
plexities  that  give  rise  to  excessive  computation  time. 

Consider  the  calculation  of  (3.10)  for  volume  scattering  using  (2.38)  in  two  di¬ 
mensions  where  r  =  (x,z): 

Zmn=  f  [  wm(r)\p'yK(r,)g0{r,r’)fn(r,)  + 

JV  Jv  L  (3.18) 

7»MV'9„(r,r')-V7„(r')]  dr  dr’ . 

One  characteristic  of  (3.18)  that  can  be  exploited  is  the  translational  invariance  of 
the  half-space  Green’s  function  g0  in  the  x-direction.  This  will  be  used  to  reduce  the 
number  of  times  go  must  be  calculated.  It  would  also  be  desirable  to  simplify  the 
integration  in  (3.18)  such  that  the  integration  of  the  Green’s  function  over  the  volume 
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is  decoupled  from  the  medium.  This  would  allow  the  computation  for  the  Green’s 
function  (for  a  specific  geometry)  to  be  performed  once  and  then  used  repeatedly 
for  different  realizations  of  a  random  medium.  These  techniques  will  be  apparent  in 
the  details  of  the  derivations  to  follow. 

In  this  section,  two  types  of  basis  and  testing  functions  will  be  used  to  approx¬ 
imate  volume  scattering.  First,  subsectional  basis  functions  that  approximate  the 
field  and  medium  as  piecewise-constant  over  a  grid  of  rectangular  cells  will  be  used. 
Then  Fourier  basis  and  testing  functions  will  be  used  as  an  example  of  the  use  of 
full-domain  functions.  The  subsectional  basis  functions  are  non-zero  only  over  a 
subsection  of  the  domain  of  T,  and  usually  provide  analytically  simpler  solutions 
than  full-domain  functions.  Alternatively,  full-domain  basis  functions  exist  over  the 
entire  domain  of  the  operator. 

With  knowledge  of  the  solution  of  a  particular  problem,  the  choice  of  full-domain 
basis  functions  can  provide  solutions  that  have  fewer  unknowns  and  are  computation¬ 
ally  more  efficient.  The  trade-off  is  usually  the  generality  of  the  numerical  method. 
In  this  chapter,  the  subsectional  basis  functions  will  provide  the  most  useful  method 
of  solving  the  half-space  volume  scattering  problem  at  hand.  All  numerical  results 
shown  will  be  found  using  this  method.  The  full-domain  solution  will  be  presented 
only  as  an  example  of  other  alternative  and  interesting  choices  that  can  reduce  the 
complexity  of  the  problem  for  specific  geometries. 

3.2.1  Point  Matching  and  Pulse  Basis  Functions 

The  most  direct  solution  using  MoM  is  found  by  choosing  testing  and  basis  functions 
that  provide  the  simplest  evaluation  of  the  scattering  matrix  elements.  One  choice 
of  testing  functions  is  to  require  that  the  solution  be  satisfied  at  discrete  points  in 
the  region  of  interest.  This  is  accomplished  by  using  Dirac  delta  functions, 


wm(x,z)  =  S(x  -  xm)5(z  -  Zm)  ; 


(3.19) 
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where  xm  and  zm  are  the  discrete  points  of  interest.  In  doing  so,  one  integral  required 
to  evaluate  the  matrix  element  is  eliminated.  This  method  is  referred  to  as  point 
matching. 

The  simplest  choice  of  basis  functions  (that  also  satisfy  orthogonality)  is  the 
pulse  basis  function, 


/»(*» z)  = 


1 

< 

0 

w 


if  |x  —  xn\  <  Ax/2  and  \z  —  zn\  <  Az/2  , 
otherwise  , 


(3.20) 


which  exists  only  over  a  small  rectangle  with  sides  Ax  and  Az  centered  at  the  nth 
discretized  point.  This  is  equivalent  to  dividing  the  scattering  volume  into  small 
cells  over  which  the  field  and  the  medium  fluctuations  are  assumed  to  be  constant: 

OO 

p{x,z)  =  J ^p{xn,zn)fn(x,z )  , 

n—  1 
oo 

z)  =  ^  ]7it(^»)^i>)/»(^)^)  >  (3.21) 

n=l 

oo 

7p(xiz)  =  ^J7P{xn,zn)fn{x,z)  * 

71=1 

The  coefficient  pn  in  (3.6)  is  simply  the  value  of  the  unknown  field  over  the  nth  cell, 

Pn=p{Xn,Zn). 

To  ensure  a  meaningful  numerical  solution,  it  is  necessary  to  employ  a  combi¬ 
nation  of  basis  and  testing  functions  that  have  finite  and  well-defined  coefficients 
throughout  the  domain  of  the  solution  [57].  This  heuristic  constraint  imposes  a  dif¬ 
ferentiability  requirement  on  the  chosen  set  of  functions  -  the  polynomial  degree  of  a 
subsectional  basis  function  must  increase  in  proportion  to  the  number  of  derivatives 
in  the  operator.  Because  the  volume  scattering  operator  contains  a  first  derivative, 
the  basis  and  testing  functions  must  have  continuous  first  derivatives. 

To  satisfy  this  requirement  and  still  be  able  to  apply  point  matching/pulse  basis 
functions,  a  weak  form  of  the  volume  scattering  operator  must  be  used.  The  first 
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derivative  of  the  field  in  (3.18)  cannot  be  eliminated,  but  can  be  approximated  with 
finite-difference  operators 


and 


p(xn,Zn )  =  2^(S*  ~  Sxl)P(Xn,Zn) 

=  ^  [P(*n+1 ,  Z)  -  p(xn.  l ,  Z)] 


-^p(Xn,Zn)  =  2^(S*  ~  S2l)p{Xn,Zn) 


(3.22) 


1 

=  2/\z  I \p(X1  Zn+l)  ~~  p{x?  5 


(3.23) 


where  sx  and  sz  and  their  inverses  are  shift  operators.  The  matrix  elements  can 
then  be  expressed  as 

Zmn  =  J  j  dx  dz  z  )9o{xmi  x  ?  Z  )  T 

cell 

7 p(x',  *')  [^°(®rn,  zm;  Xf,  z') 2^(S*  -  Sx1)  +  ^3'24^ 

~f9o(xm 5  zm',x,z) 2^ ( sz  ~  sz  )j  )  > 


and 


8 mn  —  8mn  > 


(3.25) 


where  8mn  is  the  Kronecker  delta  function.  The  vector  p0  contains  the  unperturbed 
field  points,  and  the  solution  vector  p  contains  the  unknown  field  points:  po  = 
[po {xm,zm)\  and  p  =  [p(xm,zm)\.  The  integral  in  Zmn  is  now  evaluated  over  a  cell 
centered  at  the  points  of  the  discretized  field. 

The  elements  of  the  Z  matrix  represent  the  interaction  between  cells  of  the 
discretized  field.  For  cells  far  apart  from  each  other  (off-diagonal  elements,  m  ^  n) 
it  is  expected  that  the  interaction  will  be  weak  due  to  spreading  loss  and  attenuation 
in  the  medium.  Elements  closer  to  the  matrix  diagonal  represent  the  interaction 
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between  neighboring  points  and  are  expected  to  be  larger  in  magnitude  than  the 
off-diagonal  elements. 

Rearranging  the  order  of  the  finite  difference  operators  and  redefining  the  Green’s 
function,  the  matrix  elements  take  the  form 

=  k  ^n)9oi^mi  Zmi  %ni  %n) 


7p(®w*  zn)  r\  'Eni  %n) 

uzn 

The  modified  Green’s  function  (go)  is  the  integration  of  the  half-space  Green’s  func¬ 
tion  over  the  aperture  of  the  cell: 


7 p{xni  zn)  Qq(x  m,  zm 5  xn?  zn) 


(3.26) 


QoiXmi  Zm'i  %n)  —  J ^  9o{'^mi  Zm,  &  y  Z  )  dx  dz 

cell 

1‘Xm+^T'  fZm  +  ^jr 

=  I  I  9o{xm 5  zm\  x  ,  z  )  dx  dz  . 

JXm-^L  JZm-%L 


(3.27) 


The  finite  difference  operators  in  (3.26)  operate  on  the  index  n  of  the  quantities  in 
the  square  brackets  only. 

The  integration  over  each  cell  in  the  modified  Green’s  function  (3.27)  can  be 
performed  analytically  or  numerically.  For  matrix  elements  “far  from”  the  diagonal 
elements,  the  Green’s  function  in  the  integrand  will  vary  slowly  over  each  small  cell, 
and  the  integrand  can  be  assumed  to  be  approximately  constant.  For  matrix  ele¬ 
ments  “near  to”  the  diagonal  or  on  the  diagonal,  the  integrand  of  (3.27)  cannot  be 
assumed  to  be  constant.  In  this  thesis,  a  simple  trapezoid  rule  numerical  integra¬ 
tion  was  found  to  be  sufficiently  accurate  to  approximate  the  cell  integration.  The 
definition  of  “near  to”  and  “far  from”  the  diagonal  is  a  function  of  the  behavior  of 
the  half-space  Green’s  function  and  the  cell  size.  In  general,  this  must  be  defined 
for  the  specific  geometry  and  medium  properties  of  the  problem  at  hand. 

For  diagonal  elements  (m  =  n),  the  Green’s  function  is  singular  at  the  center 
of  the  cell.  A  useful  approximation  to  solve  the  cell  integration  is  found  by  first 
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dividing  the  Green’s  function  into  the  direct  path  and  reflected  path  components, 

5„(r,r')  =  —~Ho\k(r  -  r'))  +  gr{r,r') ,  (3.28) 

where  is  the  zeroth-order  Hankel  function  of  the  first  kind.  The  reflected  path 
Green’s  function  is 

poo  1 

gr( r,  r')  =  -M  ^R(kx)e^z+zVk^x-x,)dkx  .  (3.29) ' 

J—oo  p 

The  singularity  only  exists  in  the  direct  path.  A  simple  analytic  approximation  for 
the  direct  path  is  to  replace  the  integral  of  the  direct  path  Green’s  function  over  the 
rectangular  cell  with  an  integral  over  an  equivalent  circular  cell  [29].  The  equivalent 
circular  cell  will  have  radius  a  =  y/AxAz/n.  Then,  using  the  identities 

J  rHQ1\kr)  dr  =  ^H^^kr)  (3.30) 

and 


H^’ikr)  ps  -i- 


(3.31) 


as  r  ->•  0,  where  H[^  is  the  Hankel  function  of  the  first  kind  and  first  order,  the 
diagonal  matrix  elements  can  be  approximated  as 


Az(S*  Sz  )m  T p{.xmt  zm)  9r{xmizm)  ^ 


(3.32) 


where 


(3.33) 


The  integral  over  the  cell  of  the  ^-derivative  of  the  reflected  path  and  the  x-  and  z- 
derivatives  of  the  direct  path  Green’s  functions  are  zero  due  to  symmetry.  Again,  the 
cell  integration  in  the  reflected  path  Green’s  function  can  be  performed  numerically 
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or  analytically.  However,  now  the  integration  over  the  singularity  in  the  direct  path 
is  performed  analytically. 

For  a  periodic  random  medium,  the  matrix  equations  and  the  calculation  of  the 
matrix  elements  remain  the  same,  except  for  the  Green’s  function.  The  wavenumber 
integration  in  the  half-space  Green’s  function  for  the  non-periodic  case  is  replaced 
by  a  summation  over  Floquet  modes.  The  cell  integration  is  performed  in  the  same 
manner,  and  equations  (3.26)  and  (3.32)  are  unchanged.  However,  equation  (3.27) 
is  now  defined  in  terms  of  the  periodic  Green’s  function 


)(®m?  Zm'i  ®n>  %n)  —  ^  9o (®m j  Zm'i  ®  ,  Z  )  dx  dz  . 


(3.34) 


Also,  equation  (3.33)  is  replaced  by 


Qrij^m  j  Zm )  —  J ^ j  Qr (^m 5  % m  j  ?  %  )  dx  dz  , 


(3.35) 


where  (3.29)  is  replaced  by 


oo  i 


(3.36) 


3.2.2  Fourier  Testing  and  Basis  Functions 

When  using  the  point  matching  method  to  evaluate  realistic  volume  scattering  prob¬ 
lems  the  number  of  unknowns  can  become  very  large.  An  alternative  choice  of  testing 
and  basis  functions  is  one  that  reduces  the  size  of  the  problem,  the  general  strat¬ 
egy  being  to  sacrifice  the  simplicity  of  evaluating  the  matrix  elements  in  return  for 
reducing  the  dimension  of  the  matrix  equation.  One  way  this  can  be  accomplished 
is  by  using  full-domain  basis  functions.  If  the  basis  functions  are  chosen  such  that 
the  unknown  field  is  well  defined  by  a  smaller  set  of  unknowns  (compared  to  point 
matching),  then  the  dimension  of  the  problem  can  be  greatly  reduced. 
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A  convenient  choice,  at  least  in  the  case  of  a  periodic  random  medium,  is  to  use 
Fourier  basis  and  testing  functions, 


Nk/ 2-1 

p(x,z)=  E 
l=-Nk/2 


(3.37) 


tOm(*,»)  =  WiOOe-4**-  , 

where  Afc  is  the  number  of  terms  in  the  Fourier  series.  The  field  is  expanded  in  the 
^-direction  only  where  kxi  =  2irl/ L  +  k{x,  and  Pi  is  the  depth-dependent  Fourier 
coefficient.  The  testing  function  wm  is  also  an  ^-direction  Fourier  component  with 
depth  dependence.  The  indexing  is  defined  over  two  dimensions  such  that  m  —  (i,j). 
In  addition  to  expanding  the  unknown  field,  the  medium  is  also  expanded  as  a  Fourier 


series: 


Nk/ 2-1  x 

1 «(*,*)  =  E  Mz)e,2’"I,L  , 

(3.38) 

Nk/ 2-1  V  ' 

?,(*.*)=  E 

v=-Nk/2 

For  a  periodic  medium  with  continuous  fluctuation  and  no  sharp  discontinuities,  the 
series  can  be  limited  to  a  finite  number  of  terms  (JV*)  to  sufficiently  represent  the 
unknown  field. 

Using  the  above  expansions  in  the  periodic  integral  equation  (2.116),  and  ex¬ 
ploiting  the  orthogonality  of  the  Fourier  components,  the  inner  product  integration 
of  the  method  of  moments  reduces  to  the  form 

-  aO  yO  Nk/2-l  r 


r  t ®  cO  *  r 

wm(r)Lp(r)dr  =  L2  Wi(z)dz  ^  k2A{j-i)(z')Pi(z')gj(z,z' 

Jy  J-oo  J-oo  l-_Nk/2*- 

(kxikxjPi(z')9j(z,z')  -  ~^~7 ~  — g7~)| dz'  » 


(3.39) 


where 


&(*»*0  = 


J2(*ssi)e^l*+*'1  +  e*'1*-*' 


(3.40) 


62 


In  general,  W{,  Pi,  Av,  and  Bv  are  arbitrary  functions  of  the  depth  coordinate  0.  The 
simplest  form  they  can  take,  since  the  medium  is  not  periodic  in  depth,  is  to  require 
point  matching, 

Wi(z)  =  6(z-zi)  ,  (3.41) 


and  pulse  basis  functions  in  the  ^-direction, 

Nz 

pi(z )  =  , 

k= 1 
Nz 

Av(z )  =  y]aUvfu(z)  ,  (3.42) 

U~\  •  f 

Nz 

=  buvfu(,%)  5 

where 

{1  if  | z  —  zn |  <  Az/2  , 

(3.43) 

0  otherwise. 

This  is  equivalent  to  slicing  the  medium  at  each  value  of  z  and  performing  a  Fourier 
transform  in  the  x-direction  on  each  slice. 

The  unknown  field  and  the  random  medium  are  now  expressed  as  a  two-dimensional 


grid  of  coefficients.  To  write  the  system  as  a  matrix  equation,  it  is  necessary  to  ar¬ 
range  the  unknowns  as  a  vector.  This  is  done  by  numbering  the  two-dimensional 
grid  of  coefficients  with  a  single  one-dimensional  index  defined  as  m  =  i  +  jNz  and 
n  =  k  +  INg.  Using  the  above  testing  and  basis  functions  and  the  vector  indexing 


of  the  coefficients,  the  elements  of  the  K  matrix  are  found  to  be 

7;i7i  =  I?  /\z  k  dk(j—i)9j^iizk)  kxikxjbk(j-i)9j(zi,  Zk) 
1  (a  „-l\  (h  1 

“oaJ5*  sz  )*  J  • 


(3.44) 


Again,  the  derivatives  in  the  ^-direction  have  been  approximated  using  the  finite  dif¬ 


ference  operator,  and  the  integration  over  each  cell  in  the  ^-direction  approximated 


as  in  Section  3.2.1. 
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3.3  Monte-Carlo  Simulations 

Monte-Carlo  simulations  of  volume  scattering  from  a  two  dimensional  random  medium 
are  performed  to  develop  first  and  second  moment  statistics  of  the  scattered  field. 
The  translational  invariance  of  the  Green’s  function  in  the  horizontal  direction  can 
be  exploited  to  reduce  computation  and  memory  usage.  Only  a  subset  of  the  en¬ 
tire  set  of  Green’s  function  calculations  is  performed.  Since  the  field  points  do  not 
change  position  from  realization  to  realization,  the  Green’s  function  calculations  are 
only  performed  once  for  a  set  of  realizations.  This  greatly  simplifies  the  numerical 
solution.  In  addition,  since  the  unknown  field  points  are  discretized  on  a  uniform 
grid,  the  wavenumber  integration  of  the  half-space  Green’s  function  (2.33)  can  be 
performed  efficiently  using  an  FFT  (see  Appendix  D). 


3.3.1  Incident  Field 

To  avoid  the  interaction  of  volume  elements  at  the  side  boundaries  of  a  non-periodic 
scattering  volume,  a  Gaussian  tapered  incident  field  is  used.  The  incident  field  is 
defined  at  the  interface  of  the  half-space  as  a  Gaussian  aperture  function: 


Po(x,z  =  0)  =  e~x2/92 eik°xcosei  . 


(3.45) 


The  field  below  the  interface  is  found  by  taking  the  inverse  Fourier  transform  of 
the  spectral  representation  of  the  aperture  function  and  propagating  the  field  in  the 
negative  ^-direction: 


rco 

p0(x,  z<  0)  =  -^=  /  T(kx)  e-(**-*ocosS,)V/4  eikxx-ipzdkx  # 
2v7r  J- oo 


(3.46) 


The  incident  power,  Si,  is  found  by  integrating  the  normal  component  of  the 
energy  flux  over  the  interface: 

dPo(x,z  =  0+)' 


i  r 
2poU  J- c 


Im 


p0(x,z  =  0+)- 


dz 


dx  . 


(3.47) 
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Due  to  the  orthogonality  of  Fourier  components,  the  integral  over  x  is  eliminated  and 
the  incident  power  is  found  to  be  a  function  of  a  single  integral  over  the  propagating 
wavenumbers  in  the  upper  medium: 

Si  =  I**  0o e-^-koCOse^g2/2dkx  .  (3.48) 

4^/>o  J-ko 

For  a  periodic  random  medium,  the  incident  field  will  be  taken  to  be  a  simple 
unit  amplitude  plane  wave.  The  plane  wave  field  is  incident  to  the  interface  at  the 
angle  0,-,  and  it  will  refract  into  the  volume  such  that 

po(x,z<0)=T(k„)eik-1-^  ,  (3.49) 

where  k{x  =  k0  cos  0,-,  /%  =  y/k%  —  kfx,  and  T  is  the  plane  wave  transmission  coeffi¬ 
cient.  The  incident  intensity  immediately  above  the  interface  is 

'-SSI-  (3'50) 

Figures  3.1a  and  3.2a  illustrate  typical  geometries  and  incident  fields  used  for  the 
numerical  simulations. 

3.3.2  Realization  of  a  Random  Medium 

Realizations  of  a  two-dimensional  Gaussian  random  medium  are  generated  by  Fourier 
synthesis  (see  Appendix  E).  The  assumption  of  Gaussian  statistics  is  made  only  for 
convenience.  The  medium  has  an  exponential  correlation  function  and  power-law 
spectrum,  as  described  in  Section  2.1.1.  The  correlation  function  is 

C( r)  =  a2e-r/lc  ,  (3.51) 

where  lc  is  the  correlation  length  and  cr2  is  the  variance  of  the  media  fluctuations 
(either  density  or  compressibility  fluctuations).  The  corresponding  two-dimensional 
power  law  spectrum  is 

<72/(27r/c) 

(1  HI  +  A:2)3/2  ‘ 


(3.52) 
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Typical  realizations  of  density  fluctuation  using  a  power  law  spectrum  are  shown  in 
Figures  3.1b  and  3.2b.  The  size  of  the  volume  is  constrained  by  the  incident  field 
and  statistical  convergence  of  the  numerical  solution. 

3.3.3  Statistical  Convergence  and  Accuracy 

Statistical  convergence  with  respect  to  the  random  parameters  of  the  model  must 
be  demonstrated  for  the  Monte-Carlo  simulations.  For  volume  scattering  simula¬ 
tions,  convergence  must  be  shown  with  respect  to  the  size  of  the  volume  (in  number 
of  wavelengths),  the  resolution  of  the  volume  discretization  (number  of  points  per 
wavelength),  and  the  number  of  realizations  used  to  form  ensemble  averages.  The 
minimum  size  of  the  volume  and  the  minimum  number  of  discretized  points  is  a 
function  of  the  correlation  length  of  the  medium.  A  smaller  correlation  length  may 
permit  a  smaller  volume  to  adequately  represent  the  statistics  of  the  medium;  how¬ 
ever,  it  may  require  more  discretized  points  per  wavelength.  Convergence  with 
respect  to  these  parameters  must  also  be  shown  for  a  range  of  correlation  lengths. 

Three  convergence  tests  are  performed  using  two  values  of  correlation  lengths 
(/c/ X  =  0.5  and  lc/X  =  1.0),  chosen  to  represent  typical  correlation  lengths  of  ma¬ 
rine  sediments  for  high-frequency  scattering.  Figure  3.3  shows  convergence  of  the 
numerically  calculated  total  scattering  strength  as  a  function  of  the  number  of  points 
per  wavelength  and  the  number  of  wavelengths  in  the  scattering  volume.  The  to¬ 
tal  scattering  cross-section  is  found  using  a  two-dimensional  version  of  (2.92).  The 
numerical  simulations  converge  for  the  medium  discretized  with  at  least  10  points 
per  wavelength  and  greater.  With  respect  to  the  size  of  the  volume  (Lx/ A),  the 
numerical  estimates  show  good  convergence  even  with  small  volumes:  All  of  the 
horizontal  dimensions  tested  (4  <  Lx/\  <  25)  show  good  convergence.  Figure  3.4 
shows  convergence  with  respect  to  the  number  of  realizations.  A  minimum  of  50 
realizations  is  necessary  to  form  ensemble  averages  for  the  types  of  media  used  in 
this  test  ( NxX/Lx  =  10  and  Lx/X  =  6.4). 


Incident  Field 


Figure  3.1:  Simulation  of  volume  scattering  in  a  non-periodic  medium:  (a)  The 
tapered  incident  field  \p0(x,z)[,  (b)  a  realization  of  a  Gaussian  random  medium 
7P(x,  z)\  and  (c)  the  solution  for  the  total  scattered  field  \p{x,  z) |.  The  incident  field 
is  above  critical  at  =  35°,  lc  =  0.5A,  ap  =  0.10,  \x  =  —1.1,  p/po  =  2.0,  c/cq  =  1.1, 
8  =  0.02,  with  10  points  per  wavelength.  All  coordinates  are  normalized  by  the 
wavelength. 
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Figure  3.2:  Simulation  of  volume  scattering  in  a  periodic  medium:  (a)  The  plane 
incident  field  \po(x,  z) |;  (b)  a  realization  of  a  Gaussian  random  medium  7P(x,  z);  and 
(c)  the  solution  for  the  total  scattered  field  |p(a:,  z)  |.  The  incident  field  is  subcritical 
at  6i  =  20°,  /c  =  1.0A,  o’ p  =  0.10,  p  =  -2.0,  p/p0  =  2.0,  c/c0  =  1.1,  8  =  0.02,  with 
10  points  per  wavelength.  All  coordinates  are  normalized  by  the  wavelength. 
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The  accuracy  of  the  simulation  method  can  be  evaluated  by  comparing  the  nu¬ 
merical  results  with  the  solution  of  problems  that  can  be  solved  analytically.  The 
comparison  will  test  applicability  of  the  choice  of  basis  functions,  the  accuracy  of  the 
numerical  evaluation  of  the  half-space  Green’s  function,  and  the  accuracy  of  the  cell 
integration  methods.  Simulating  the  phase  matching  conditions  required  to  model 
perfect  reflection  and  transmission  in  a  non-random  medium  will  provide  an  upper 
bound  on  the  numerical  errors  of  the  random  medium  simulations. 

Consider  the  transmission  of  an  incident  plane  wave  into  an  infinite  half-space 
with  a  linear  profile  (in  depth)  of  the  square  of  the  index  of  refraction  (7*  =  az)  and 
constant  density  profile  (7 p  =  b).  The  wave  equation  for  this  case  can  be  written  as 

V2$(x,z)  +  P  $(*,*)  =  0,  (3.53) 

where 

«(*,*)=  (3.54) 

yPo 

The  solution  for  the  field  in  the  LHS  is  well  known  [11,  49],  and  is  expressed  as  an 
Airy  function 

$(x,z)  =  W(kx)Ai[t{z)]eik*x  ,  ^  (3.55) 


where 


t(z )  =  —  2/#, 
H  =  (ak2)~1/3, 


W(kx) 


t0  =  H2[k2x-k2] 

_ 2ip(30H 

pQhi'(to)  +  ip(3ahi{to)H 


(3.56) 


Figure  3.5  shows  the  numerical  results  of  the  test  case  for  a  non-periodic  medium 
with  a  normally  incident  tapered  field.  The  analytic  solution  for  the  tapered  field  is 
found  by  performing  a  plane  wave  decomposition  and  applying  (3.55)  to  each  Fourier 
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component  of  the  field.  Figure  3.6  shows  the  results  of  the  test  case  for  a  periodic 
medium  with  a  plane  incident  wave.  The  sound  speed  gradient  in  both  cases  is 
defined  with  c/co  =  1.2,  p/p0  =  2.0,  a  =  —0.05  and  b  =  0.1.  The  piecewise  constant 
curves  (stepped  curves)  are  the  numerical  MoM  solutions.  The  solid  curves  are  the 
analytic  results,  and  the  dashed  curves  are  the  incident  fields.  A  discretization  of 
10  points  per  wavelength  was  used  in  the  simulations.  The  MoM  solution  compares 
well  with  the  analytic  results  for  both  cases. 

3.3.4  Matrix  Inversion 

In  this  thesis,  an  iterative  conjugate  gradient  method  was  used  to  perform  the  ma¬ 
trix  inversion  [58].  No  approximations  were  used  to  reduce  the  size  of  the  matrix. 
For  the  problems  considered,  this  standard  method  was  sufficiently  accurate  and 
computationally  efficient.  A  scattering  matrix  for  sediment  volume  scattering  is 
typically  8192  x  8192  matrix  elements  in  size  (corresponding  to  a  volume  of  length 
12. 8A  and  depth  6.4A  discretized  with  10  points  per  wavelength).  For  larger  prob¬ 
lems,  sparse  matrix  methods  or  other  numerical  approximations  that  exploit  the 
loss  in  the  medium  in  order  to  reduce  the  size  of  the  scattering  matrix  could  be 
considered  [18,  78].  However,  the  effort  required  to  improve  the  numerical  efficiency 
of  a  particular  solution  must  be  weighed  against  the  ever-increasing  computational 
resources  that  are  becoming  available  each  year  at  decreasing  costs. 

An  alternative  to  matrix  inversion  techniques  is  to  use  the  numerical  small  per¬ 
turbation  method  of  Section  3.1.1.  Some  problems  that  lie  outside  the  range  of 
validity  of  first-order  perturbation  theory  can  still  be  solved  with  sufficient  accuracy 
using  higher-order  perturbation  approximations.  The  perturbation  solution  can  be 
found  to  any  higher  order  by  efficient  matrix  multiplications  using  (3.14).  An  ex¬ 
ample  is  shown  in  Figure  3.7.  The  scattering  strength  (as  defined  in  Section  2.5)  is 
found  using  both  the  numerical  perturbation  method  and  by  matrix  inversion  (using 
the  conjugate  gradient  method).  Each  higher-order  perturbation  solution  (plotted  as 


Total  Scattering  Strength  (dB)  Total  Scattering  Strength  (dB) 
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Number  of  Points  per  Wavelength  (NXA/LX) 


Number  of  Wavelengths  (Lx/A) 


Figure  3.3:  Monte-Carlo  convergence  with  respect  to:  (a)  the  number  of  points  per 
wavelength;  and  (b)  the  size  of  the  medium.  The  symbol  *  represents  lc  =  0.2A,  the 
o  represents  lc  =  0.5A,  and  the  x  represents  lc  =  1.0A.  A  total  of  50  realizations 
was  used  to  estimate  the  ensemble  average. 


Total  Scattering  Strength  (dB) 


Figure  3.5:  Fields  within  a  half-space  medium  with  a  linear  wavenumber  profile 
(c/co  =  1.2,  p/po  =  2.0,  a  =  —0.05,  b  =  0.1)  for  a  tapered  incident  field:  (a) 
Contour  plot  of  magnitude  of  the  tapered  field  MoM  solution;  (b)  Real  component  of 
the  MoM  solution  (stepped  curve)  compared  with  the  analytic  solution  (solid  curve) 
and  the  incident  field  (dashed  curve);  (c)  Imaginary  components  of  the  solutions. 
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(a)  (b) 


Re[p(x=0,z)]  lm[p(x=0,z)] 


Figure  3.6:  Fields  within  a  half-space  medium  with  a  linear  wavenumber  profile 
(c/co  =  1.2,  p/p0  =  2.0,  a  =  —0.05,  b  =  0.1)  for  a  periodic  medium  and  plane  wave 
incident  field:  (a)  Real  components  of  the  MoM  solution  (stepped  curve)  compared 
with  the  analytic  solution  (solid  curve)  and  the  incident  field  (dashed  curve);  (b) 
Imaginary  components  of  the  solutions. 
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the  solid  dots)  is  compared  with  the  exact  solution  (plotted  as  the  circles).  Note  that 
since  the  random  medium  is  Gaussian  the  odd-order  contributions  to  the  scattering 
strength  are  zero.  The  perturbation  series  must  be  taken  to  fourth-order  to  model 
the  scattering  strength  accurately.  An  ensemble  of  200  realizations  of  the  periodic 
random  medium  was  used  to  find  the  scattering  strength.  If  a  problem  lies  outside 
the  range  of  validity  of  the  perturbation  method  altogether,  the  perturbation  series 
may  not  converge  towards  the  correct  solution.  Higher-order  perturbation  methods 
must  be  used  with  caution. 

3.4  Summary  and  Conclusions 

Two-dimensional  acoustic  volume  scattering  from  a  random  medium  can  be  simu¬ 
lated  numerically  using  the  Method  of  Moments.  Compressional  wave  attenuation 
in  the  medium  limits  the  size  of  the  scattering  volume.  Therefore,  the  solution  do¬ 
main  is  limited  in  size,  and  numerical  solutions  can  be  performed  without  the  use 
of  excessive  computational  resources  and  time. 

Two  MoM  solution  methods  are  presented,  in  which  the  medium  is  modeled 
as  both  periodic  and  non-periodic.  A  tapered  incident  field  is  used  with  the  non¬ 
periodic  medium  in  order  to  limit  the  interaction  of  the  incident  field  with  the 
edges  of  the  volume.  A  plane-wave  incident  field  is  used  with  the  periodic  medium. 
Periodicity  is  assumed  only  in  the  horizontal  direction.  In  both  cases,  the  vertical 
extent  (depth)  of  the  volume  is  defined  by  the  distance  the  incident  field  penetrates 
into  the  lossy  medium. 

Point  matching  and  pulse  basis  functions  are  demonstrated  for  use  with  both 
the  periodic  and  non-periodic  volumes.  Full-domain  Fourier  basis  and  testing  func¬ 
tions  can  be  used  with  a  periodic  medium.  In  both  cases,  the  scattering  matrix 
is  diagonally  dominant  and  potentially  sparse.  Convergence  and  accuracy  tests  for 
Monte-Carlo  simulations  show  that  the  medium  should  be  minimally  discretized  at 
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(a)  (b) 


0  50  100  150  0  50  100  150 
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Figure  3.7:  Comparisons  of  MoM  solutions  for  the  volume  scattering  strength  found 
using  the  iterative  SPM  method  and  the  exact  matrix  inversion  method.  The  circles 
(o)  are  the  exact  solutions,  and  the  solid  dots  (•)  are  the  SPM  solutions  including 
terms  up  the  specified  order:  (a)  second-order  SPM  using  (2.82);  (b)  third-order 
SPM  using  (2.84);  (c)  fourth-order  SPM  using  (2.86);  and  (4)  fifth-order  SPM.  The 
incident  field  is  subcritical  at  0;  =  10°,  lc  —  1.0A,  ap  =  0.10,  p  =  —2.0,  p/ po  =  2.0, 
c/co  =  1.1,  8  =  0.02,  with  10  points  per  wavelength. 
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10  points  per  wavelength  in  order  to  model  scattering  due  to  the  medium  fluctua¬ 
tions.  The  minimal  size  of  the  scattering  volume  required  for  convergence  is  smaller 
than  expected.  Monte-Carlo  convergence  was  demonstrated  for  a  horizontal  dimen¬ 
sion  as  small  as  Lx/\  =  4.  Random  realizations  of  the  medium  are  generated  by 
Fourier  transforms  assuming  Gaussian  statistics.  The  number  of  realizations  re¬ 
quired  to  estimate  the  ensemble  average  of  the  scattered  field  should  be  >  50  to 
assure  statistical  convergence. 
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Chapter  4 

VALIDITY  OF  THE  SMALL  PERTURBATION  METHOD 

Insight  into  determining  whether  single  or  multiple  scattering  is  significant  in  a 
continuous  random  medium  can  be  gained  by  examining  the  scales  and  interrelations 
of  various  parameters  of  the  problem.  The  most  essential  of  these  parameters  are 
the  wavelength  of  the  incident  field,  the  variance  of  the  fluctuations  in  the  medium, 
the  correlation  length  of  the  fluctuations,  and  the  characteristic  length  scale  over 
which  scattering  occurs  in  the  volume  (due  to  attenuation  in  the  medium). 

In  general,  it  is  desirable  to  have  analytic  relations  and  criteria  for  predicting  the 
occurrence  of  multiple  scattering.  For  the  general  acoustic  fluid  half-space  problem, 
however,  this  is  difficult  and  not  at  present  available.  Criteria  would  have  to  prop¬ 
erly  account  for  the  interface  and  reflection  back  into  the  medium  of  the  scattered 
field.  The  effects  of  correlation  of  the  medium  compressibility  and  density  fluctua¬ 
tions  must  also  be  considered.  An  alternative  analysis  can  be  made  using  numerical 
simulations.  By  comparing  solutions  of  a  weak  scattering  method,  such  as  per¬ 
turbation  theory,  with  exact  numerical  methods  that  do  not  make  weak  scattering 
assumptions,  the  range  of  validity  of  the  perturbation  method  can  be  inferred. 

In  this  chapter,  both  multiple  scattering  theory  and  numerical  methods  (as  in 
Chapter  3.1)  will  be  used  to  investigate  the  applicability  of  perturbation  methods  for 
volume  scattering  from  a  fluid  sediment  half-space.  Parameters  of  the  random  media 
are  limited  to  those  of  marine  sand  and  soft-mud  sediments  modeled  as  fluids.  The 
physical  parameters  of  the  sediment  (e.g.,  mean  density,  sound  speed,  and  statistical 
properties  of  the  medium)  are  selected  from  current  publications  of  studies  made  of 
sediment  core  samples  and  in  situ  measurements  [60,  23,  13]. 
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4.1  Multiple  Scattering  Theory 

Classical  multiple  scattering  theory  involves  using  diagram  techniques  to  sum  the 
perturbative  series  for  the  mean  field  (Dyson’s  equations)  and  the  intensity  or  covari¬ 
ance  of  the  field  (Bethe-Salpeter  equation)  in  the  scattering  medium.  The  Bethe- 
Salpeter  equation  is  directly  related  to  the  radiative  transfer  equation  where  attenu¬ 
ation  in  the  medium  is  a  function  of  the  intrinsic  absorption  of  the  medium  and  the 
scattering  attenuation  [80].  Scattering  attenuation  is  due  to  incident  energy  diverted 
from  the  incident  direction  by  multiple  scattering. 

Radiative  transfer  theory  provides  an  analytic  method  of  solving  for  the  effects 
of  multiple  scattering  on  the  field  intensity,  where,  in  general,  there  is  no  solution  to 
the  Bethe-Salpeter  equation.  Radiative  transfer,  on  the  other  hand,  is  a  phenomeno¬ 
logical  approach  that  does  not  provide  insight  into  the  scattering  processes  and  the 
parametric  relationships  of  the  medium  and  the  incident  wave  that  dictate  the  oc¬ 
currence  of  multiple  scattering.  To  gain  physical  insight  into  multiple  scattering  it 
is  more  useful  to  study  Dyson’s  equation.  To  proceed  analytically,  assumptions  are 
typically  made,  the  most  basic  of  which  is  the  bilocal  approximation. 

A  great  deal  of  literature  is  available  on  multiple  scattering  theory.  The  bulk  of 
this  work  focuses  on  scattering  due  to  the  index  of  refraction  (created  by  permit¬ 
tivity  fluctuations)  in  an  unbounded  medium.  The  sediment  problem  differs  in  that 
the  medium  is  a  lossy  half-space.  This  constrains  the  incident  and  scattered  fields 
to  a  region  near  the  interface.  The  sediment  problem  also  differs  in  that  density 
fluctuations,  as  well  as  compressibility  fluctuations,  must  be  considered.  In  general, 
bistatic  scattering  into  the  UHS  is  of  greater  interest  than  forward  scattering  into 
the  LHS.  Therefore,  the  correlation  properties  of  the  density  and  compressibility 
may  have  a  significant  role  in  suppressing  or  enhancing  multiple  scattering  effects. 

Effective  medium  parameters  for  an  unbounded,  statistically  homogeneous  fluid 
with  strong  density  and  compressibility  fluctuations  have  been  discussed  by  Zhuck 
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[89,  90]  in  the  low-frequency  limit  ( k  0).  This  work  extends  the  weak  scattering 
criterion  developed  in  previous  multiple  scattering  theories,  which  states 

|*b  L  <t*|2  <  1  ,  (4.1) 

to  include  density  fluctuations  that  obey  the  same  relation, 

|*b  Ip  <r„|2  «  1  .  (4.2) 

The  length  scales  lK  and  lp  are  the  correlation  lengths  of  the  compressibility  and 
density  fluctuations,  respectively.  The  cross-correlation  of  density  and  compressibil¬ 
ity  fluctuations  is  not  considered.  In  addition,  larger  correlation  scales  (comparable 
to  the  acoustic  wavelength)  and  loss  due  to  intrinsic  attenuation  are  not  considered. 

The  role  of  half-space  effects  in  multiple  scattering  has  been  investigated  using 
Dyson’s  equation  to  find  a  renormalized  effective  wavenumber  of  the  mean  field 
[43,  63].  The  principal  conclusion  was  that  the  mean  field  in  the  bounded  medium 
cannot  be  modeled  by  defining  a  single  effective  wavenumber,  as  with  the  unbounded 
medium.  Instead,  a  transition  layer  near  the  boundary  must  be  defined  in  which  the 
incoherent  wave  reflected  from  the  interface  affects  the  coherent  wave  to  the  same 
order  as  an  incoherent  wave  generated  without  the  interface. 

In  this  section,  the  bilocal  approximation  will  be  applied  to  the  acoustic  scat¬ 
tering  case  where  compressibility  and  density  fluctuations  are  present  and  of  com¬ 
parable  magnitude.  The  general  acoustic  scattering  problem,  including  half-space 
effects,  will  not  be  considered.  Compressional  wave  attenuation  will  be  included  as 
well  as  the  effects  of  density  and  compressibility  correlation.  This  differs  from  devel¬ 
opments  in  the  literature  where  only  compressibility  (or  permittivity)  fluctuations 
are  considered.  In  general,  the  problem  of  scattering  due  to  density  (or  permeability) 
fluctuations  has  not  been  widely  studied.  In  electromagnetics,  strong  permeability 
fluctuations  are  rarely  encountered.  In  ocean  acoustics,  multiple  scattering  is  usu¬ 
ally  only  considered  in  the  forward  direction  where  sound  speed  fluctuations  are  the 
dominant  scattering  mechanism  and  density  fluctuations  need  not  be  considered. 
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(4.3) 


4-1.1  Dyson’s  Equation 

Consider  the  small  perturbation  solution  for  the  acoustic  field  from  Section  2.3 
rewritten  as 

CO 

p( r)  =  g( r,  r0)  =  <to(r,  r0)  +  Y  9n(r,  r0)  , 

n= 1 

where  g  is  the  field  at  r  due  to  a  source  at  r0,  and 

gn+i( r,r0)  =  J  t)go(r, r')gn(r', r0)  + 

7p(r') VVo(r,  r')  •  Vgn (r',  r0)] dr'  . 


(4.4) 


Assume  (for  convenience  only,  since  measurements  neither  support  nor  challenge  this 
assumption)  that  the  fluctuations  in  the  medium  are  zero-mean  Gaussian  random 
variables  such  that  the  even  moments  are  given  by 


(7>c(rl)  ■  ’  *  7«(r2n))  —  ^  ^  CKK(ri,  Fj)  ■  •  •  CKK{ Ffc,  F/) 

pairs 

(7«(ri)  ■  •  •  lp{r2n))  =  Y  C*p(ri’  ri)  * ' '  C*p(rk,  r/) 

pairs 

(7t(ri)-"7e(*2n))  =  ^2  C„(ri,rj)--C„(n,ri) 


(4.5) 


pairs 


(4.6) 


and  the  odd  moments  are  zero.  For  example, 

(7«(ri)7#>(r2)7K(r3)7p(r4))  =CKp(ri,  r2)C'«/)(r3,  r4)  +  C'K/t(r1,r3)C'pp(r2,r4) 

+  C'Kp(ri,r4)C'pK(r2,r3)  . 

The  renormalized  perturbation  series  for  the  mean  field  (g)  is  formed  by  performing 
an  ensemble  average  over  realizations  of  the  medium: 

g(r,r0)  =  (g(r,r0)) 

OO 

=  5'o(r,r0)-f  ^(5-„(r,r0)) 

71=1 

=  9o(r,  r,j)  +  (fc(r,  rQ))  +  (gt(r,  r0))  +  •  • ;  . 


(4.7) 
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Odd-order  terms  are  zero  due  to  the  Gaussian  assumption. 

Solving  this  system  of  equations  by  successive  iterations,  Dyson’s  equation  for  a 
random  medium  with  both  density  and  compressibility  fluctuations  is  found  to  be 


g{  r,  r0)  =<7o(r,ro)  +  k4  /M  ,r1)mKK(r1,r2)^(r2,ro)dridr2 
+  k2  JJ £o(r,  ri)mKp(ri,  r2)  •  V2<j(r2,  r0)drxdr2 
+*JJ  Vi5f0(r,  rx)  •  mpK(ru r2)g(r2, r0)drxdr2 
+  JJ  Viflf0(r,ri)  .mPp(ri,r2)  •  V2£(r2,  r0)cMr2  . 


(4.8) 


The  operator  Vx  is  the  gradient  with  respect  to  rx;  V2  is  the  gradient  with  respect 
to  r2;  and  V2VX  is  a  dyad  operator.  This  form  of  Dyson’s  equation  is  similar  to  the 
standard  Dyson’s  equation  found  in  [75,  68],  except  the  mass  operator  m( rx,r2)  now 
ha§  four  components  that  represent  the  interactions  of  scattering  from  density  and 
compressibility  fluctuations.  Also,  note  that  the  integration  limits  in  (4.8)  are  over 
the  half-space  volume.  For  now,  the  effects  of  the  interface  will  be  ignored  and  the 
limits  of  integration  will  be  extended  over  infinity.  Consequently,  the  unperturbed 
free-space  Green’s  function  can  be  used  for  go  instead  of  (2.33). 

The  mass  operators  are  best  described  using  the  diagram  method  introduced  in 
Section  2.3.2  with  the  addition  of  several  new  symbols.  Represent  the  correlations 
between  density  and  compressibility  fluctuation  with  a  dashed  line: 


k4CKK(ri,Tj)  ~  • - •  , 

*i  rj 

k  CKp(Vi)  Tj)  ^  •  O  ? 

r*  Tj 

k2CpK( ritrj)  ~  o - •  , 

r*  r  j 

Cpp(Ti,  Tj)  ~  O - O  . 


(4.9) 


Represent  the  mean  Green’s  function  for  the  inhomogeneous  medium  with  a  double 


82 


line: 


§(r,  p0) 


r  r0 


(4.10) 


The  mass  operators  are  represented  as  double  dashed  line  segments: 

mKK(ri,ri)  ~  •===•  , 

*1  r2 

m«p(ri,r2)  ~  •===°  , 

ri  r2 

m„(ri,r2)  ~  o===»  , 


(4.11) 


ri  r2 


mpp{ ri,r2)  ~  o===o  . 

ri  V2 

Following  the  procedure  of  separating  weak  and  strong  scattering  interactions  [68], 
Dyson’s  equation  for  an  acoustic  fluid  can  then  be  written  in  diagram  representation 
as 

=  +  - 
=:  H — 


+ 

+ 


where  the  mass  operators  are  defined  as  follows: 


/  /  N  N  /  /  v 

>===•  =  i- - ^  +  » - » - • - •  +  • - * - « - • 


(4.12) 


#===o  =  m- - £>  +  • - o - • - o  +  • - o - •  c 
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Consider  a  statistically  homogeneous  medium  where  the  correlation  functions, 
the  Green’s  functions,  and  the  mass  operators  are  all  functions  of  difference  coordi- 
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nates  only.  Dyson’s  equation  is  then  rewritten  as  a  double  convolution  integral: 

g{ r  -  r0)  =go(r  -  r0)  +  k*  JJ  go{r  -  ri)mKfC(ri  -  r2)p(r2  -  r0)dridr2 
+  P  JJ tfo(r  -  r1)mKp(r1  -  r2)  •  V2^(r2  -  r0)dridr2 
+  P  JJ  Vifif0(r  -  ri)  •  mpK(ri  -  r2)^(r2  -  r0)dr1dr2 

+!!  Vi#0(r  -  rx)  •  mpp( ri  -  r2)  •  V2flf(r2  -  r0)dridr2  . 

The  unperturbed  free-space  Green’s  function  in  three  dimensions  is 


(4.13) 


9o(r)  = 


Jkr 


47r  r 


(4.14) 


Solving  for  the  mean  Green’s  function  using  Fourier  transforms  with  the  transform 
of  the  unperturbed  Green’s  functions, 

1 


8tt3(A;2  -  P)  ’ 

and  the  transform  of  the  mass  operators, 


Me 


l(i(k)  =  J  rn(rl  -  r2)„(je  ,k''V/r  , 


(4.15) 


(4.16) 


an  expression  for  g  can  be  found  explicitly  as 


§(r-r°)  =  8 


eik-(r-r0) 


-dk 


k2  -  M( k) 

The  function  M  is  found  in  terms  of  the  Fourier  transforms  of  each  component, 


(4.17) 


M(k)  =  iSM.Jk)  -  2iPk  •  M.,(k)  -  k  •  M„(k)  ■  k  ,  (4.18) 

A 

with  MKp  =  MpK.  The  function  MKK  is  a  scalar,  MKp  is  a  vector  pointing  in  the  k 

A  A 

direction,  and  Mpp  is  a  dyad  whose  orientation  is  defined  by  kk. 

Equation  (4.17)  is  the  mean  Green’s  function  for  the  coherent  field.  In  this  form, 
(4.17)  is  exact.  The  expression  for  the  mass  operator,  however,  is  an  infinite  expan¬ 
sion  and  is  generally  not  obtainable  in  closed  form.  Dyson’s  equation  is  therefore 
only  useful  in  that  it  allows  approximations  to  be  made  of  the  mass  operator  that 
lead  to  insight  into  the  behavior  of  g. 
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4-1.2  Bilocal  Approximation 

The  simplest  approximation  to  Dyson’s  equation  is  the  bilocal  approximation  which 
involves  keeping  only  the  first  terms  in  the  mass  operators  [10],  such  that 

mK«(ri,ri)  »  k4CKK(rur2)go(rur2)  , 
m«p(ri,r2)  «  P^(ri,r2)V2So(ri,r2)  , 

(4.19) 

mpK(r1,r2)  «  fc2^K(r1,r2)Vit/0(ri,r2)  , 
mpp(YUT2)  »  C'pp(ri,r2)ViV25fo(ri5r2)  • 

The  bilocal  approximation  is  represented  in  diagram  notation  as 

=  =  -  +  - 1 - V=  +  - l - *= 

+  — <^=  + _ 


(4.20) 


Applying  this  approximation  to  (4.8),  the  mean  field  is 
g(r,r0)  =go(r, r0)  +  k4  j J CKK(rur2)g0(r, ri)g0(r1,r2)g(r2,r0)dr1dr2 

k2  /I  CKp(ru  r2)g0(r,  r1)V2fif0(ri,  r2)  •  V2#(r2,  r0)cMr2 

k2  II  CpK(ri,r2)Vig0(r,r1)  •  Vi^0(ri,r2)^(r2,r0)dr1dr2 

I j  Cpp(ri,r2)Viflb(r,ri)  •  ViV2flb(ri,r2)  •  V2^(r2,r0)dr1dr2  . 

Again,  consider  a  statistically  homogeneous  medium  where  the  density  and  com¬ 
pressibility  fluctuations  are  proportional  using  (2.13)  and  (2.14).  Using  the  bilocal 
approximation,  the  individual  mass  operators  that  comprise  (4.18)  can  be  written 


+  . 
+ 
+ 


as 


Mkk( k)  »  g2 1  Cpp(r)g0(r)e~ik'rdr  , 

M*p(k)  «  g  I  Cpp(r)Vg0(r)e~lk'rdr  ,  (4.21) 

Mpp  (k)  «  I  Cpp(r)VVgo(r)e-ik-rdv  . 
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Using  the  free  space  Green’s  function  in  three  dimensions  with  an  isotropic  and 
exponential  correlation  function, 

C„( r)  =  ,  (4.22) 

each  term  in  the  mass  operator  can  be  derived  analytically  as 

<r2/2' 

M«„(k)  =  ,  (4.23) 


k  •  M„p(k)  =  ,kl<  (c  -  _  (2i>  -  ac) 


and 

4-^)1 

? 

where 

a  =  1  —  iklc  , 

6  =  iklc  , 

c  =  ln[(a  +  b)/(a  —  6)]  . 


cr2k2 

k-M„(k)-k=-^r 


ikL 


(2b  —  ac)  —  (a  - 


+  (26  —  ac)(a  —  b2) 


(4.24) 


(4.25) 


(4.26) 


Since  equations  (4.23),  (4.24),  and  (4.25)  are  not  found  in  previous  literature,  the 
details  of  performing  the  integrations  are  outlined  in  Appendix  F. 

For  an  isotropic  medium,  (4.17)  can  be  rewritten  in  spherical  coordinates  with 
vector  k  pointing  in  the  ^-direction.  Performing  the  integration  over  the  angle 
variables,  a  single  integral  results: 


If  eikR 

9 ^  =  4; fHR  J  F-£2-M(fc)fc</A:  ’ 

— OO 


(4.27) 
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with  R  =  |r  — r0|.  Taking  residues  of  the  integrand  with  respect  to  k  at  the  dominant 
pole,  the  mean  Green’s  function  for  the  medium  can  be  found  approximately  as 


m 


eikeR 

4nR 


(4.28) 


It  is  expected  that  a  single  pole  ( ke  «  k)  in  the  upper  half  of  the  complex  k  plane 
will  dominate  the  total  solution  of  (4.27).  The  other  wave  solutions  -  corresponding 
to  the  other  poles  -  are  expected  to  decay  rapidly  and  can  be  ignored.  This  has 
been  shown  for  the  case  of  index  of  refraction  fluctuation  only  [68],  and  is  assumed 
to  be  true  for  the  current  case  if  the  medium  fluctuations  are  not  too  strong  and 
|M(*)|<£. 

The  wavenumber  ke  is  the  effective  wavenumber  of  the  renormalized  medium. 
The  real  component  of  ke  yields  the  effective  phase  speed  of  the  medium.  As  the 
mean  path  of  the  wave  through  the  medium  is  increased  by  multiple  scattering, 
the  real  part  of  ke  increases,  thus  causing  an  effective  decrease  in  the  phase  speed. 
The  imaginary  part  of  ke  corresponds  to  the  total  coherent  field  attenuation  in  the 
effective  medium.  Multiple  scatter  will  act  to  increase  the  total  attenuation  due  to 
scattering  loss  out  of  the  incident  direction. 

The  onset  of  multiple  scattering  can  be  examined  by  studying  the  behavior  of 
the  dominant  pole  ke  as  the  parameters  of  the  medium  are  varied.  The  poles  of 
(4.27)  are  the  roots  of  the  denominator  of  the  integrand: 


k2-k2-  M(k)  =  0  .  (4.29) 

In  the  weak  scattering  limit,  the  position  of  ke  should  approach  that  of  Jc  in  the 
complex  plane.  Therefore,  rather  than  solving  for  the  root  exactly,  using  (4.23), 
(4.24)  and  (4.25),  the  position  of  the  dominant  pole  ke  can  be  found  approximately 
by  taking  a2  to  be  small  and  using  perturbation  analysis.  In  the  zeroth-order  solution 
for  the  pole,  where  there  is  no  multiple  scattering,  M{k)  =  0  and  ke  =  k.  The  first- 
order  solution  is  then  found  by  substituting  the  zeroth-order  result  into  (4.29),  such 
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that 


k2e  -P-  M(k )  «  0  . 


(4.30) 


With  M(k)  small,  the  effective  wavenumber  for  the  the  bilocal  approximation  is 


ke 


«  k  + 


M(k ) 
2k  ’ 


(4.31) 


Multiple  scattering  can  be  examined  by  plotting  the  changes  in  the  real  and 
imaginary  components  of  ke,  defined  respectively  as 


,  Ref&e  —  fc] 
e  =  Re[k] 


(4.32) 


and 


Im[fce  —  k] 

Im[fc] 


(4.33) 


Figures  4.1  and  4.2  show  contours  of  7'  and  7"  as  functions  of  lc  and  fi  for  several 
values  of  the  attenuation  coefficient  ol\.  All  plots  are  made  with  a 2  =  0.01.  The 
dashed  contours  in  Figure  4.2  are  values  of  7''  found  by  solving  for  the  roots  of 
(4.29)  numerically  using  a  complex  simplex  method.  The  solid  contours  are  values 
of  7"  found  using  ke  estimated  with  (4.31).  The  perturbation  method  for  finding  the 
effective  wavenumber  shows  reasonable  correspondence  with  the  numerical  method 
for  small  lc/X  and  fi  near  —1.  When  using  the  perturbation  method,  7'  and  7"  are 
(to  first  order)  proportional  to  the  variance  of  the  medium.  Therefore,  it  is  only 
necessary  to  plot  contours  for  a  single  value  of  a2.  Contours  for  other  values  of  a2 
will  scale  proportionally. 

The  range  of  parameters  shown  in  Figures  4.1  and  4.2  reflects  typical  values 
observed  in  marine  sediment  over  a  variety  of  environments  and  sediment  types. 
The  plot  for  a\  =  0.1  (which  corresponds  to  typical  mud  and  silt  sediment)  shows 
higher  scattering  attenuation  of  the  mean  field.  As  the  intrinsic  attenuation  in  the 
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medium  (Im[&])  increases,  multiple  scattering  decreases  and  the  scattering  atten¬ 
uation  (Im[fce  —  k])  decreases.  A  value  of  ct\  =  1.0  corresponds  to  a  typical  sand 
sediment. 

Relative  to  the  homogeneous  wavenumber,  changes  in  the  real  part  of  the  effective 
wavenumber  are  small  as  expected  (7'  <C  7").  Changes  in  the  imaginary  part  of  the 
effective  wavenumber  are  a  minimum  for  values  of  /z  ~  —1.  This  can  be  attributed 
to  the  minimum  in  the  sound  speed  fluctuations  when  density  and  compressibility 
fluctuations  are  negatively  correlated. 

4.1.3  First-order  Multiple  Scattering 

For  weak  scattering  -  but  not  necessarily  single  scattering  -  the  imaginary  component 
of  ke  that  is  due  to  scattering  attenuation  should  approach  the  total  scattering 
cross-section  derived  from  the  first-order  perturbation  approximation.  The  effective 
wavenumber  found  using  (4.31)  can  then  be  compared  to  the  first-order  perturbation 
result  of  Section  2.2.  In  other  words,  the  power  attenuation  of  the  mean  field  should 
obey  the  relation 

2Im[A:e  —  k]  =  atv  =  /  crv  du  ,  (4.34) 

J  47T 

where  av  is  the  differential  volume  scattering  cross-section  and  du  =  sin  Qd0d<f).  This 
approach  is  often  referred  to  as  first-order  multiple  scattering  [30]. 

From  Section  2.3,  using  (2.90)  and  neglecting  half-space  effects,  the  first-order 
differential  cross-section  is 

ff«  =  |l*IV  +  o»  WSM,  (4.35) 

where  0  is  the  angle  between  the  incidence  and  scattering  directions,  and  kd  = 
2Re[fe]  sin(#/2)  is  the  Bragg  wavenumber.  Using  the  correlation  function  (4.22),  the 
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total  volume  scattering  cross-section  is 


<*tv  -  cr]ll\k\ 


2p2  2p> 

(  2 ab  \ 

a2  -  b2  b2  1 

^  a2  —b2  J 

+ 


K26- 


2ac  + 


2a2b 
a 2  —  b2 


(4.36) 


where  now  a,  b,  and  c  are  redefined  as 


a  =  1  +  2(Re[fc]fc)2  , 


b  =  — 2(Re[fc]/c)2  , 


c  =  In 


(4.37) 


Figure  4.3  shows  contour  plots  of  the  ratio  7"  using  ke  found  from  (4.36)  and  (4.34) 
for  various  sediment  types  (solid  contours).  Also  plotted  in  Figures  4.3  (as  the 
nearly  hidden  dashed  contours)  is  7"  found  using  (4.31).  A  small  difference  in  the 
two  methods  of  calculating  ke  can  be  seen  at  lc/ A  1  near  p  =  —  1. 


4.2  A  Comparison  with  Numerical  Methods 

Multiple  scattering  theory  provides  an  analytic  method  of  studying  higher-order 
scattering  in  an  unbounded  medium.  For  a  random  half-space,  however,  analytic 
results  are  not  easily  obtained.  The  half-space  geometry  prohibits  a  solution  by 
Fourier  transforms  because  the  medium  and  the  Green’s  function  are  no  longer 
translationally  invariant  in  the  z-direction.  The  integration  over  the  medium  is  no 
longer  an  integration  over  an  infinite  continuum,  but  rather  over  the  half-space. 
Analytical  results  cannot  be  found  unless  further  approximations  are  made. 

Numerical  methods  will  be  employed  to  analyze  half-space  effects  and  their  effect 
on  multiple  scattering.  Two  potential  effects  will  be  investigated:  1)  the  existence 
of  a  transition  layer  near  the  interface  as  predicted  by  [43,  63]  and  its  effect  on 
the  definition  of  an  effective  wavenumber  in  the  sediment;  and  2)  the  effect  of  the 
half-space  and  multiple  scattering  on  the  volume  scattering  cross-section. 
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Figure  4.3:  First-order  multiple  scattering  compared  to  the  bilocal  approximation: 
Contours  of  7^  as  a  function  of  /C/A,  //,  and  cv,\ .  The  solid  contour  lines  are  found 
using  first-order  multiple  scattering,  the  dashed  contours  (nearly  hidden  below  the 
solid  contours)  are  found  using  the  bilocal  approximation  with  a  small  variance.  All 
contours  are  calculated  with  crp  —  0.01. 
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Simulations  are  performed  in  two  dimensions,  and  the  multiple  scattering  theory 
is  for  three  dimensions.  Therefore,  the  comparison  of  the  results  of  numerical  simula¬ 
tions  with  the  results  of  multiple  scattering  theory  is  limited  to  relative  magnitudes. 
The  exact  magnitudes  of  the  effective  wavenumbers  found  using  both  methods  can 
not  be  compared  directly.  Multiple  scattering  theory  could  be  derived  for  two- 
dimensional  scattering.  However,  the  derivation  of  the  mass  operator  is  much  more 
difficult  with  the  two-dimensional  Green’s  function  (a  Hankel  function),  and  it  is 
not  done  in  this  thesis. 

4-2.1  Half-Space  Effects  on  the  Mean  Field 

Consider  first  the  scattered  field  within  the  sediment.  As  in  the  previous  section 
where  Dyson’s  equation  was  used,  multiple  scattering  can  be  observed  by  examining 
the  mean  field.  The  mean  field  is  generated  by  performing  Monte-Carlo  simulations 
and  averaging  the  total  field  over  a  sufficient  number  of  realizations  to  remove  the 
incoherent  field.  A  periodic  medium  is  used  to  allow  for  plane  wave  propagation,  and 
the  incident  field  is  normal  to  the  interface.  Six  cases  are  examined  using  parameters 
of  the  medium  that  correspond  to  typical  sand  sediments.  The  cases  are  defined  in 
Figure  4.4  over  a  range  of  p  and  correlation  scales,  lcf  A.  Multiple  scattering  theory 
predicts  that  in  these  cases  the  effective  wavenumber  will  be  influenced  by  at  least 
first-order  multiple  scattering.  Figure  4.5  illustrates  the  geometry,  the  incident  field, 
and  mean  fields  within  the  scattering  volume  of  the  numerical  simulations.  The 
incident  field  (Figure  4.5a)  is  normal  to  the  interface  at  z  =  0.  The  mean  field 
(Figure  4.5b)  is  an  average  over  50  realizations. 

The  bilocal  approximation  for  an  unbounded  medium  predicts  an  effective  wavenum¬ 
ber  at  which  the  mean  field  propagates  in  the  scattering  medium.  However,  near 
the  interface  the  definition  of  an  effective  medium  is  obscured.  Figure  4.6  shows 
numerically  estimated  values  of  7"  as  a  function  of  depth  into  the  half-space,  found 
using  the  estimated  effective  wave  number  of  the  mean  field.  The  estimates  of  the 
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Figure  4.5:  The  magnitude  of:  (a)  the  incident  field;  (b)  the  mean  field;  and  (c) 
the  horizontally  averaged  real  components  of  the  incident  field  (dashed  line)  and 
mean  field  (solid  line)  within  the  scattering  volume  of  a  periodic  medium.  Sediment 
parameters  correspond  to  case  B  in  Figure  4.4.  The  incident  field  is  normal  to  the 
interface  at  z  =  0.  All  coordinates  are  normalized  by  the  wavelength. 
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effective  wavenumber  are  found  by  solving  for  the  ^-component  of  the  exponent  of 
the  wave  field  at  each  point  in  the  depth  direction: 

[ke  -  k]  •  £  =  ln^V_P°l  .  (4.38) 

—iz 

The  estimate  of  the  mean  field  (p)  is  found  by  averaging  Monte-Carlo  simulations, 
and  po  is  the  incident  field.  Within  a  transition  layer  near  the  interface,  the  scattering 
attenuation  increases  from  approximately  zero  at  the  interface  to  a  constant  value. 
The  depth  at  which  the  effective  wavenumber  becomes  constant  is  the  depth  of 
the  transition  layer.  In  all  cases,  the  depth  of  the  transition  layer  is  larger  than 
the  correlation  length  of  the  medium,  with  larger  correlation  lengths  showing  larger 
transition  layers. 

As  mentioned  previously,  no  direct  comparison  can  be  made  between  the  numer¬ 
ically  estimated  effective  wavenumber  and  the  three-dimensional  bilocal  approxima¬ 
tion.  However,  the  relative  magnitudes  of  the  test  cases  agree  with  the  theoretical 
predictions. 

J.2.2  Bistatic  Scattering  Strength 

Next,  consider  the  scattered  far-held  in  the  UHS  by  estimating  the  equivalent  surface 
scattering  strength  for  the  medium.  Again,  a  periodic  medium  is  assumed.  The 
scattering  strength  is  found  by  estimating  the  variance  of  the  scattered  far-held  for 
a  large  number  of  realizations.  Because  a  periodic  medium  is  used,  the  scattering 
angles  are  limited  by  the  Floquet  modes,  and  the  scattering  cross-section  is  dehned 
using  (2.133).  Figure  4.7  shows  the  bistatic  scattering  strength  for  different  incident 
angles  for  case  B  of  Figure  4.4.  A  critical  angle  exists  at  24.6° 

The  bistatic  scattering  strengths  for  various  incident  angles  above  and  below  the 
critical  angle  show  good  agreement  with  hrst-order  perturbation  results.  The  solid 
and  dashed  lines  are  the  formally  averaged  perturbation  theory  with  and  without 
half-space  effects.  The  dots  are  the  numerically  calculated  scattering  strength.  Half- 
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space  effects  are  significant  in  the  specular  direction.  Although  multiple  scattering 
is  observed  in  the  scattering  medium  (as  shown  by  the  mean  field),  the  scattering 
strength  is  modeled  well  by  first-order  perturbation  methods. 

4-2.3  Assessment  of  Multiple  Scattering  in  Marine  Sediments 

Figure  4.8  shows  four  different  contour  plots  of  7''  for  four  different  marine  sites. 
The  sediment  physical  properties  at  each  site  were  estimated  from  cores  samples 
[34,  14].  The  measured  sediment  parameters  p  and  Zc/A  are  plotted  as  points  (*)  on 
the  contour  plots.  For  the  estimated  sediment  parameters  and  the  frequencies  used, 
multiple  scattering  theory  predicts  that  no  significant  multiple  scattering  of  the  field 
within  the  sediment  will  occur  at  three  of  the  sites  (7^  <  5%).  The  Panama  City 
site  shows  fairly  strong  scattering  within  the  sediment  (7"  10%).  The  sediments 

at  the  Eckernfoerde  and  Orcas  sites  are  primarily  silt  and  mud.  The  Key  West  and 
Panama  City  sites  are  primarily  sand. 

A  primary  indicator  of  multiple  scattering  in  marine  sediments  (at  least  at  these 
four  sites)  is  the  amplitude  of  the  cross-correlation  between  the  compressibility  and 
density.  In  general,  compressibility  and  density  are  negatively  correlated  and  of 
comparable  magnitude.  This  leads  to  p  1,  and  multiple  scattering  is  minimal. 
For  values  of  p  away  from  negative  unity,  compressibility  and  density  act  to  en¬ 
hance  multiple  scattering.  The  Panama  City  site  in  Figure  4.8  shows  non-negligible 
multiple  scattering  due  to  the  relatively  high  estimate  of  p  &  —2.5. 

4.3  Summary  and  Conclusions 

The  bilocal, approximation  to  Dyson’s  equation  for  the  coherent  field  is  extended  to 
include  density,  as  well  as  compressibility,  fluctuations  in  volume  scattering  from  a 
continuous  Gaussian  random  medium.  The  Gaussian  medium  assumption  is  only 
used  to  develop  a  bilocal  approximation  for  sediment  volume  scattering  that  can  be 
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used  to  study  the  applicability  of  the  small  perturbation  method.  It  is  not  intended 
as  a  realistic  model  of  sediment  fluctuations. 

The  effective  wavenumber  in  the  scattering  medium  is  found  using  the  bilocal 
approximation.  Alternatively,  the  effective  wavenumber  is  found  using  a  first-order 
multiple  scattering  method.  The  two  expressions  for  the  effective  wavenumber  agree 
well  for  typical  values  of  sediment  physical  properties.  However,  they  differ  slightly 
at  lc/\  <C  1  near  //  =  —  1  (as  shown  numerically). 

The  effective  wavenumber  within  the  sediment  half-space  is  found  by  numerical 
Monte-Carlo  simulations  using  the  method  of  moments  (for  several  test  cases).  The 
numerical  simulations  show  the  existence  of  a  transition  layer  near  the  sediment- 
water  interface,  as  predicted  by  one-dimensional  theory  [63,  43].  The  UHS  scattering 
strength  for  a  typical  marine  sediment  is  also  calculated  using  numerical  simulations. 
The  bilocal  approximation  for  the  effective  wavenumber  in  the  sediment  predicts 
that  multiple  scattering  is  present  in  some  cases.  However,  numerical  simulations 
of  the  scattering  strength  in  the  UHS  still  show  good  agreement  with  the  first-order 
perturbation  model  (2.88)  when  half-space  effects  are  included. 

One  conclusion  that  can  be  drawn  from  this  analysis  is  that  good  estimates  of 
the  physical  properties  of  sediment  are  critical  for  making  predictions  of  scattering. 
Good  estimates  of  fx  are  especially  important.  The  current  measurements  of  fi 
indicate  that  first-order  multiple  scattering  within  the  sediment  is  significant  for  the 
Panama  City  site,  for  example.  However,  the  core  analysis  used  in  this  case  may 
not  have  provided  adequately  accurate  estimates  of  (jl  to  make  reliable  predictions. 

A  final  judgment  as  to  the  applicability  of  the  weak  scattering  assumption  in 
marine  sediments  cannot  be  made.  Further  research  is  needed  to  address  the  pro¬ 
portionality  of  sediment  density  and  compressibility  fluctuations.  If  the  assumption 
of  proportionality  can  be  verified,  the  range  of  values  of  ji  for  typical  sediments 
must  be  established.  If  a  constitutive  relation  between  density  and  compressibility 
cannot  be  established,  the  effects  of  uncorrelated  scattering  between  density  and 
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compressibility  fluctuations  must  be  addressed. 

For  far-field  scattering  from  the  sediment  into  the  UHS,  an  important  conclusion 
can  be  made.  In  general,  it  is  observed  from  numerical  simulations  that  even  when 
multiple  scattering  is  present  in  the  medium,  first-order  estimates  for  the  scattering 
strength  are  good  approximations  to  the  true  scattering  strength.  This  is  assumed 
to  be  the  case  in  typical  marine  sediments,  and  it  is  assumed  to  be  true  for  the  cases 
shown  in  Figure  4.8. 

An  explanation  for  this  is  most  likely  found  in  the  transition  layer  near  the 
sediment-water  interface.  As  observed  through  numerical  simulations,  multiple  scat¬ 
tering  diminishes  near  the  interface.  This  is  shown  in  the  estimates  of  the  depth- 
dependent  effective  wavenumber  in  the  scattering  medium  (Figure  4.6).  The  scat¬ 
tered  field  in  the  UHS,  however,  is  more  strongly  influenced  by  scattering  from  the 
sediment  near  the  interface  -  because  attenuation  limits  the  depth  to  which  the 
incident  field  will  penetrate.  Therefore,  if  multiple  scattering  is  diminished  in  the 
volume  of  the  sediment  near  the  interface  (which  most  influences  the  scattered  far- 
field),  it  is  expected  that  weak  scattering  methods  will  model  the  scattered  far-field 
adequately.  In  other  words,  half-space  effects  that  create  a  transition  layer  in  the 
upper  sediment  diminish  the  effects  of  multiple  scattering  on  the  far-field  scattering 
strength. 
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Chapter  5 

BIOLOGICAL  SOURCES  OF  HETEROGENEITY  IN 

SEDIMENTS 

Identifying  the  specific  mechanisms  by  which  benthic  organisms  mix  sediment 
is  critical  in  the  development  of  a  mathematical  model  of  bioturbation.  However, 
generalization  to  an  arbitrary  level  of  description  seems  exceedingly  difficult  due  to 
the  variety,  complexity  and  irregularity  of  bioturbational  processes.  Exceptions  to 
the  general  behavior  of  organisms  will  always  require  unique  treatment  in  a  model. 
For  example,  communities  of  organisms  may  affect  actions  of  individuals  in  ways 
that  may  be  unique  to  a  particular  environment  [4].  It  seems  realistic,  however,  to 
model  the  cumulative  effects  of  bioturbation  and  the  subsequent  effects  on  sediment 
properties  in  a  statistical  sense  -  that  is,  to  develop  models  that  describe  the  stochas¬ 
tic  properties  of  the  mixing  process  without  attempting  to  model  the  behavior  of 
individual  organisms. 

This  chapter  will  begin  with  a  brief  review  of  the  activities  of  benthic  organisms 
that  cause  sediment  mixing  and  heterogeneity  (at  least  the  activities  that  affect  the 
physical  properties  of  sediments).  Then,  a  review  of  current  mathematical  models 
of  bioturbation  is  presented,  keeping  in  mind  that  these  models  were  developed  for 
other  purposes  (and  other  parameterizations)  than  for  acoustic  remote  sensing.  Fi¬ 
nally,  a  new  model  is  presented  that  closely  follows  the  existing  models,  but  which  is 
extended  to  include  higher-order  statistical  measures  of  fluctuations  in  the  sediment 
bulk  properties.  The  fluctuation  statistics  can  then  be  used  to  relate  the  parameters 
of  the  mixing  model  to  the  statistics  of  the  scattered  sound  field. 
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5*1  Mixing  of  Sediments  by  Benthic  Organisms 

The  activities  of  organisms  (meiofauna  and  macrofauna)  that  live  and  feed  in  marine 
sediments  rework  the  sediment  and  cause  mixing.  These  activities  are  characterized 
as  deposit  feeding,  locomotion ,  and  home  building  [83]  which  in  turn  are  characterized 
as  either  causing  local  or  nonlocal  displacement  of  sediment  particles  and  pore  water 
[1,  9].  Local  displacement  is  a  continuum  of  small  random  movements  that  transport 
sediment  short  distances.  Nonlocal  displacements  are  large  (relative  to  animal  body 
size)  and  discrete  movements  of  sediment.  In  either  case,  displacements  can  take 
place  in  the  horizontal  and  vertical  directions  with  time  and  distance  scales  that 
may  vary  with  direction  and  depth  into  the  sediment.  The  precise  meaning  and 
scale  of  local  and  nonlocal  will  be  defined  later. 

Deposit  feeding  is  the  bulk  ingestion  and  subsequent  egestion  of  particles  within 
or  upon  the  sediment.  It  is  the  dominant  type  of  feeding  activity  that  displaces 
sediment  particles.  In  general,  most  deposit  feeders  must  ingest  large  amounts  of 
sediment  material  -  several  times  their  body  weight  per  day.  They  often  move 
the  sediment  distances  comparable  to,  or  greater  than,  the  length  of  the  animal. 
Depending  on  the  body  size  of  the  animal  in  question  and  the  mixing  length  scales  of 
interest,  deposit  feeding  may  produce  either  local  or  nonlocal  mixing.  Some  deposit 
feeders  restrict  their  activities  to  either  horizontal  or  vertical  mixing,  potentially 
creating  anisotropic  mixing.  Conveyer  belt  or  tunnel  feeders,  which  make  up  the 
majority  of  deposit  feeders,  tend  to  move  sediment  from  some  depth  in  the  sediment 
to  the  surface.  Other  animals  transport  sediment  in  nearly  equal  distances  in  both 
directions. 

Locomotion  is  the  movement  of  an  animal  from  one  location  to  another  within 
the  volume  or  on  the  surface  of  the  sediment  which  results  in  the  sediment  being 
displaced  in  some  manner.  Movements  may  be  due  to  feeding  or  predator-prey 
interactions,  for  example.  Since  a  large  amount  of  energy  must  be  expended  during 
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movements  within  the  volume  of  the  sediment,  infaunal  mixing  of  this  nature  will 
have  relatively  small  particle  displacements  (local  mixing)  and  be  less  intensive  than 
mixing  by  deposit  feeding.  Epifaunal  locomotion,  on  the  other  hand,  is  less  difficult, 
and  mobile  epifauna  are  generally  more  abundant.  This  potentially  creates  higher 
localized  horizontal  mixing  rates  at  the  sediment  interface.  It  may  also  influence 
vertical  mixing  by  the  “piping”  of  surface  sediments  into  open  tubes  and  burrows 
that  extend  vertically  into  the  sediment. 

Home  building  is  the  creation  of  tubes  and  burrows  for  dwellings  in  the  sediment. 
The  size  and  structure  of  tubes  and  burrows  reflect  the  size  of  the  animals  inhabiting 
them  and  the  environment  in  which  they  are  built.  Large  volumes  of  sediment  can 
be  displaced  and  moved  nonlocally  by  home  building  -  comparable  to  mixing  caused 
by  deposit  feeding.  Mixing  is  further  intensified  as  uninhabited  homes  collapse  or 
become  filled  as  sediment  falls  into  them. 

Other  factors  that  affect  mixing  (that  cannot  be  simply  classified  as  one  of  the 
above)  are  lumped  into  a  general  category  of  mixing  termed  incidental  movements. 
These  include  mixing  that  cannot  be  directly  associated  with  an  animal  activity  but 
which  occurs  as  an  indirect  result  of  the  above  activities.  For  example,  through  a 
siphoning  action  a  surface  deposit-feeder  will  frequently  dislodge  sediment  particles 
which  it  does  not  ingest  [83]. 

5.2  Classical  Diagenetic  Models 

Diagenesis  refers  to  the  processes  that  bring  about  changes  in  time  and  space  of  a 
sediment  or  sedimentary  rock  after  deposition.  Bioturbation  is  an  important  aspect 
of  this  process  in  marine  sediments  and  has  received  a  significant  amount  of  attention 
and  modeling.  The  current  models  typically  describe  the  net  effects  of  bioturbation 
as  localized  diffusion  or  biodiffusion.  This  section  is  a  brief  overview  of  the  most 
widely  published  model  [4,  9,  83]  and  provides  a  starting  point  for  the  development 
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of  a  more  general  mixing  model. 

Quantitative  analysis  of  diagenesis  requires  the  use  of  mathematical  models  that 
describe  the  mixing  process.  Classical  diagenetic  models  describe  the  steady-state 
vertical  mixing  of  the  mean  concentration  of  a  sediment  tracer  (usually  a  chemical 
or  radioisotope  tracer)  as  an  idealized  mass  balance  between  adjacent  layers  of  sedi¬ 
ment.  Using  the  sediment-water  interface  as  the  origin,  the  rate  of  change  of  a  solid 
or  liquid  tracer  at  a  fixed  depth  is  expressed  as 


+  R  > 


(5.1) 


dC_  _  _dF_ 

dt  dz 

where  C  is  the  bulk  concentration  of  a  biodiffusing  tracer  (to  be  identified).  The 
units  of  C  are  mass  per  unit  volume  of  sediment.  The  function  F  is  the  flux  of  C 
in  the  vertical  direction,  with  dimensions  of  mass  per  unit  area  per  unit  time.  The 
R  term  is  the  rate  of  diagenetic  reaction  affecting  the  concentration  or  the  rate  of 
removal  or  addition  of  tracer  from  outside  the  control  volume. 

For  loose  marine  sediments,  this  type  of  rate  equation  is  used  to  describe  the 
conservation  of  fluid  and/or  solid  trace  species.  The  tracer  may  refer  to  the  concen¬ 
tration  of  a  dissolved  solute  in  the  sediment  pore  water,  such  that 


C  =  4C,t 


(5.2) 


where  Cj  is  mass  of  solute  per  unit  volume  of  pore  water,  and  4>  is  the  sediment 
porosity.  The  concentration  may  also  describe  a  solid  species,  in  which  case 

C  =  (1  -  <f>)Cs  ,  (5.3) 


where  Cs  is  the  mass  of  tracer  per  unit  volume  of  solid.  Solids  are  typically  discrete 
sediment  grains  immersed  in  the  pore  water.  The  trace  quantity  may  also  refer 
to  the  total  fluid  mass,  in  which  case  Cj  =  <f>p/ ,  or  total  solid  mass  such  that 
Cs  =  (1  -  <f>)ps,  where  p/  and  ps  are  the  fluid  and  solid  mass  densities,  respectively. 
The  trace  quantity  may  also  be  the  total  combined  sediment  mass,  such  that 


C  =  (ftp/  +  (1  -  4>)ps  • 


(5.4) 
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The  flux  F  is  the  result  of  diffusion  and/or  advection.  The  advection  of  a  tracer 
relative  to  the  sediment-water  interface  is  primarily  due  to  depositional  burying, 
compaction,  and,  depending  on  the  environment,  externally  impressed  hydrological 
flow.  Deposition  results  in  a  net  flux  downward  into  the  sediment  relative  to  the 
sediment-water  interface.  Compaction  of  sediment  due  to  gravity  results  in  a  vertical 
flow  of  pore  water.  On  the  time  scales  of  interest  in  this  thesis,  deposition  and 
compaction  will  be  neglected.  In  principle,  nonlocal  mixing  can  also  cause  advection. 
For  example,  deposit  feeding  may  result  in  a  net  upward  or  downward  advection  of 
sediment  due  to  “conveyer-belt”  feeding. 

Local  diffusional  flux  is  in  general  due  to  bioturbation  and  molecular  diffusion. 
Considering  only  bioturbational  diffusion,  and  ignoring  molecular  diffusion,  the  to¬ 
tal  localized  flux  within  the  sediment  is  proportional  to  the  gradient  of  the  tracer 
concentration.  The  net  effect  of  bioturbation  is  lumped  into  a  single  parameter,  the 
biodiffusion  coefficient  Db,  such  that 


F  = 


(5.5) 


The  diffusive  nature  of  bioturbation  is  not  the  result  of  the  activities  of  a  single  type 
of  organism,  but  it  is  rather  the  net  effect  of  a  collection  of  infauna  having  a  wide 
range  of  mixing  lengths  that  are  small  compared  to  the  length  scales  of  the  volume 
of  sediment  being  considered. 

The  resulting  vertical  biodiffusion  equation  is 


dC  d2C 
dt  Dbdz 2 


=  R  . 


(5.6) 


A  typical  model  assumes  a  zone  of  biodiffusional  mixing  near  the  sediment-water 
interface  (0  >  z  >  — L )  over  which  Db  is  constant  and  below  which  animal  activities 
are  negligible  with  Db  =  0.  The  depth  at  which  bioturbational  mixing  stops  is 
referred  to  as  the  rework  depth  L.  Below  the  rework  depth,  other  flux  mechanisms 
are  dominant,  such  as  molecular  diffusion.  To  model  anisotropic  and  inhomogeneous 
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biodiffusion,  the  diffusion  equation  can  be  generalized  into  three  dimensions, 

^C(r,i)-V-[A(r,()VC(r,()]  =  JJ(r,(),  (5.7) 

where  the  biodiffusion  coefficient  in  now  a  dyad. 


5.2.1  Nonlocal  Mixing 


Nonlocal  mixing  is  defined  as  any  mixing  that  cannot  be  described  as  diffusive.  To 
complete  this  definition  it  is  necessary  to  define  what  is  meant  by  diffusive  in  the 
context  of  sediment  mixing.  To  start,  consider  conservation  of  a  trace  quantity  ct 
at  layer  i  in  a  one-dimensional  sediment  comprised  of  thin  layers  [9,  8].  The  rate  of 
change  of  e;  due  to  the  instantaneous  exchange  of  sediment  from  all  others  layers  is 

N 

(5.8) 


^  =  Y.  K<>c>Az  -  E  Az  ■ 


dt 


3= 1 


j=l 


The  constant  specifies  the  rate  at  which  sediment  is  exchanged  from  layer  j 
to  layer  i  (with  units  of  inverse  time  and  distance).  Figure  5.1  illustrates  exchange 
between  a  finite  number  of  layers.  Deposition  at  the  sediment-water  interface  creates 
advection  downward  relative  the  the  interface-centered  coordinate  system.  Local 
exchange  is  illustrated  as  exchange  between  adjacent  layers,  and  nonlocal  exchange 
is  between  non-adjacent  layers.  Generalizing  to  a  continuum  in  three  dimensions, 
(5.8)  is  rewritten  as  an  integro-differential  equation  in  terms  of  a  general  exchange 
function  K  that  describes  all  scales  of  mixing: 


^rc(r,t)  =  f  K(r,r'-,t)c(r',t)dr' 

dt  JVr  (5.9) 

—  c(r,t)  j  K(r',r-,t)dr' . 

Jv 

Conservation  of  the  trace  quantity  follows  from  (5.9)  by  integrating  over  the  control 
volume: 


7t  lc{rJ)dr  = 


o . 


(5.10) 


Figure  5.1:  Local  and  non-local  mixing  between  layers  of  sediment. 
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In  general,  sediment  is  composed  of  a  collection  of  complex  and  heterogeneous 
sediment  grains.  In  describing  a  continuous  trace  quantity  (concentration,  density, 
etc.)  it  is  necessary  to  consider  a  sample  volume  over  which  the  point  function  (e.g 
mass  of  particular  substance)  is  spatially  averaged.  This  is  accomplished  by  defining 
the  sample  filtering  function  s(r),  such  that 

C(r,t)=  f  c(r',t)s(r  —  r')dr' ,  (5.11) 

Jv 

where 


(5.12) 


The  sampling  filter  must  have  a  characteristic  dimension  (width)  large  enough  to 
define  the  specific  quantity  of  interest  but  small  enough  to  provide  adequate  spatial 
resolution.  Usually,  the  scale  of  the  required  spatial  resolution  dictates  the  sample 
filtering  scale.  For  example,  acoustic  waves  interact  with  fluctuations  in  the  sediment 
bulk  properties  at  length  scales  comparable  to  the  acoustic  wavelength.  Therefore, 
an  appropriate  sample  filter  for  measuring  sediment  acoustic  properties  (density, 
sound  speed)  should  have  a  length  scale  much  greater  than  the  grain  size  but  smaller 
than  the  acoustic  wavelength. 

To  gain  insight  into  the  effects  of  sample  filtering  on  the  mixing  equation,  the 
operator  (5.11)  is  applied  to  the  mixing  equation  (5.9)  to  yield 


C(r,t)  =  JJ  K(r",  r';  t)c(r',  i)s(r  -  v")dr'dr" 

V 

-  JJ  K(r',r";t)c(v",t)s(r  -  i •")*'*"  . 


(5.13) 


A  comparison  of  the  spatial  scales  of  the  exchange  function  and  the  filtering  func¬ 
tion  will  give  definition  to  the  terms  local  and  nonlocal  mixing  [Jackson-personal 
communication]. 

Local  mixing  is  defined  as  exchange  on  scales  smaller  than  the  scale  of  the  filtering 
function  -  that  is,  when  K(r",  r';i)  restricts  mixing  to  length  scales  that  are  much 


Ill 


smaller  than  the  width  of  s(r  —  r").  If  this  is  assumed,  the  filtering  function  can  be 
expanded  in  a  Taylor  series  in  (r"  —  r'): 


s(r  -  r")  =s(r  -  r')  +  (r'  -  r")  •  Vs(r  -  r') 
+  j[(r'  -  r")  ■  v]2s(r  -  r')  +  •• 


(5.14) 


(5.15) 


Using  the  first  three  terms  in  the  Taylor  series,  and  assuming  K{ r"  ,r')  restricts 
|r'  —  r"|  to  scales  much  smaller  than  the  width  of  s,  then  s(r  —  r')  ~  s(r  —  r")  within 
the  mixing  integrals.  The  mixing  equation  (5.13)  can  now  be  rewritten  as 

^C(r,  t)=  Jv  u(r'i  *)  •  v  [c(r'>  *Mr  “  r')] dr>  + 

V  •  f  D(r',  t )  •  V  [e(r',  t)s(r  —  r')]  dr'  , 

Jv 

where  the  new  mixing  parameters  are  defined  as 

u(r;,  t)  =  J  K(r",  r';  t)(r"  —  r')dr" 

and 


(5.16) 


D(r',  t)  =  l  f  K{r",  r';  i)(r"  -  r')2dr"  .  (5.17) 

2  Jv 

The  function  D  is  the  dyadic  diffusion  coefficient,  and  u  is  a  vector  function  that  de¬ 
scribes  advection  in  the  local  mixing  process  due  to  asymmetry  in  the  local  exchange 
function  about  the  point  r.  If  u  and  D  are  constants  in  space,  or  slowly  varying  func¬ 
tions  compared  to  the  filtering  function,  then  (5.15)  reduces  to  a  diffusion  equation 
in  terms  of  the  filtered  quantity: 

jtC( r.  i)  =  u(r, i)  •  VC( r,  t)  +  D( r,  ()V2C(r,  i)  .  (5.18) 

If  the  exchange  function  is  an  even  function  of  space,  such  that 


K( r)  =  K(- r)  , 


(5.19) 
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then  the  coefficient  u  will  vanish,  and  there  will  be  no  advection.  Assuming  isotropy 
in  the  exchange  function,  the  familiar  diffusion  equation  results: 


J^(r,i)=D(r,i)V2C”(r,().  (5.20) 

Therefore,  local  mixing  is  defined  as  sediment  exchange  on  scales  smaller  than  the 
sample  filter  length  scale.  Local  exchange  is  characterized  by  diffusive  mixing. 

Nonlocal  mixing  is  defined  as  exchange  of  sediment  on  scales  much  larger  than 
the  sample  filter  length  scale.  Nonlocal  mixing  is  non-diffusive.  In  this  regime 
of  mixing,  the  filter  function  s(r)  is  a  narrow  function  in  space  compared  to  the 
exchange  function  AT(r).  The  mixing  equation  (5.13)  can  then  be  rewritten  in  terms 
of  the  averaged  quantity  C ,  and  the  nonlocal  mixing  equation  results: 


Jn7(r,  t)  =  J  K(r,  r';  t)C(r',  t)dr' 

—  C(r,t)  f  K(r',r;t)dr'  . 

Jv 


(5.21) 


5.3  Inhomogeneous  Biodiffusion  of  Sediment  Bulk  Properties 


The  mixing  of  a  sediment  tracer  (for  example,  a  radioisotope  with  a  relatively  long 
half-life)  is  modeled  remarkedly  well  by  the  homogeneous  and  steady-state  solution 
of  the  biodiffusion  equation.  This  has  been  shown  by  comparison  of  theoretical 
results  with  measured  vertical  profiles  of  tracer  concentration  taken  from  sediment 
cores.  Typically,  many  cores  are  taken  in  a  particular  area  to  develop  a  statistical 
ensemble.  Each  core  profile  represents  a  single  realization  of  the  random  mixing 
process.  Any  particular  core  may  include  disturbances  associated  with  sediment 
microenvironments  such  as  strong  gradients  in  porosity  around  burrows  or  corpses 
of  larger  organisms  [4,  9].  The  ensemble  of  core  profiles  is  then  averaged  to  find  the 
mean  concentration  of  the  tracer  of  interest  as  a  function  of  depth.  The  average 
profile  is  then  compared  with  the  theoretical  profile  for  a  biodiffusive  process  and  a 
biodiffusion  coefficient  is  inferred. 
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In  the  context  of  acoustic  scattering,  however,  the  fluctuations  in  the  cores  from 
realization  to  realization  are  of  primary  interest.  Acoustic  scattering  provides  a 
measure  of  the  correlation  structure  of  the  sediment  density  fluctuations.  Therefore, 
to  use  acoustic  scattering  as  a  means  of  measuring  bioturbational  mixing,  it  is  nec¬ 
essary  to  investigate  the  transient  and  heterogeneous  nature  of  the  mixing  process 
and  develop  a  model  that  can  be  used  to  describe  higher-order  statistics  -  at  least 
the  correlation  of  the  medium  fluctuations. 

A  starting  point  to  develop  such  a  model  is  the  mixing  equation  (5.21).  Since 
this  model  is  intended  for  acoustic  remote  sensing,  the  trace  species  of  interest  will 
be  the  sediment  bulk  density.  Therefore,  using  (5.4),  with 

C  =  p  =  4>Pf  +  (1  —  <t>)p*  >  (5.22) 

and  using  the  mixing  equation  (5.21),  the  rate  of  change  of  total  mass  (fluid  and 
solid)  is  described  as 

Jj/K r,  t)  =  K(r',  r;  t)p{ r',  t)dr'  -  p( r,  t)  J  K( r,  r';  t)dr'  .  (5.23) 

By  using  the  bulk  density  as  the  trace  quantity,  it  is  assumed  that  the  solid  and 
fluid  mass  of  the  sediment  grains  and  pore  water  are  mixed  as  a  single  quantity. 
Alternatively,  the  porosity  could  be  considered  constant,  and  only  the  sediment 
grains  (without  pore  water)  would  be  subject  to  biodiffusional  mixing  [4].  This 
macroscopic  view  of  mixing  ignores  the  particular  mechanism  by  which  an  organism 
transports  sediment.  The  proper  trace  quantity  to  use  for  a  particular  sediment 
type  and  environment  has  yet  to  be  established.  Only  bulk  density  will  used  in  this 
thesis;  however,  the  modeling  methodology  is  general. 

To  start,  several  assumptions  about  the  mixing  process  are  made.  First,  owing 
to  the  modal  characteristic  of  the  size  spectrum  of  organisms,  it  is  assumed  that 
mixing  takes  place  on  two  scales:  1)  local  biodiffusional  mixing  that  represents 
the  continuous  diffusive  activity  of  smaller  organisms  (meiofauna)  or  the  collapsing 
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and  filling  of  borrows  created  by  larger  organisms;  and  2)  nonlocal  mixing  that 
represents  the  activities  of  larger  organisms  (macrofauna).  On  the  time  scales  of 
interest,  the  movement  of  sediment  is  not  necessarily  restricted  to  localized  mixing 
between  adjacent  horizontal  layers.  Rather,  larger  organisms  may  transfer  sediment 
between  non-adjacent  layers  creating  nonlocal  (non-diffusive)  mixing  that  gives  rise 
to  spatial  heterogeneity.  Furthermore,  over  short  time  periods,  it  cannot  be  assumed 
that  a  large  collection  of  organisms  with  a  broad  spectrum  of  sizes  has  reworked 
the  sediment.  Nonlocal  mixing  may  be  dominated  by  a  single  or  several  types  of 
organisms  with  characteristic  length  scales.  Therefore,  the  second  assumption  is 
that  the  activities  of  the  larger  macrofaunal  organisms  are  discrete  in  both  space 
and  time,  and  the  activity  of  meiofauna  is  assumed  to  be  continuous  in  space  and 
time. 

The  exchange  function  K  describes  the  spatial  and  temporal  characteristics  of. 
mixing.  In  the  simplest  form  (as  a  continuous  function  or  a  constant),  it  describes 
total  mixing  within  the  sediment  volume.  .  In  this  case,  every  point  in  the  volume  ex¬ 
changes  sediment  with  every  other  point  in  the  volume;  mixing  occurs  on  all  scales. 
To  describe  a  two-scale  mixing  process  (as  motivated  above)  K  is  divided  in  two 
components  -  local  and  nonlocal  exchange  -  that  have  a  finite  (but  not  overlap¬ 
ping)  spatial  extent.  The  components  are  assumed  to  be  independent  of  each  other; 
therefore,  they  can  be  superimposed  to  describe  the  total  exchange  as 

K( r',  r;  t)  =  Ki( r',  r;  t )  +  Kni( r',  r;  t )  .  (5.24) 

The  function  Ki  represents  the  continuous  and  random  local  exchange  of  sediment 
by  meiofauna  (diffusive  mixing).  The  function  Kni  represents  the  discrete  (in  time 
and  space)  nonlocal  mixing  due  to  larger  animals. 

The  random  and  discrete  nature  of  the  nonlocal  mixing  is  described  by  modeling 
Kni  as  a  series  of  source  and  sink  functions  that  are  distributed  randomly  in  space 
and  time.  To  start,  consider  only  two  such  events  -  a  source  and  sink  event  related 
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Figure  5.2:  Discrete  nonlocal  mixing  between  two  layered  regions  of  sediment.  Sed¬ 
iment  is  removed  from  a  volume  defined  by  the  sink  function  h  and  deposited  in  a 
volume  defined  by  the  source  function  g. 


to  the  removal  and  deposition  of  sediment  by  a  single  animal  due  to  a  single  action 
(a  single  feeding  motion,  for  example).  Suppose  that  a  random  mass  of  sediment  q 
is  removed  (nonlocally)  at  time  £  from  a  location  described  by  the  function  h(r  —  a). 
The  position  vector  a  is  at  the  center  of  the  volume,  and  h  describes  the  shape  of  the 
volume  of  sediment  where  the  mass  was  removed.  The  sediment  is  then  deposited 
at  some  later  time,  £,  in  a  location  defined  by  the  function  g(r  —  b),  where  b  is  the 
central  location  and  g  describes  the  shape  of  the  volume  of  sediment  where  the  mass 
is  deposited.  While  the  sediment  is  being  transferred,  assume  that  it  is  temporarily 
removed  from  the  system  and  no  mass  is  lost  during  transport.  The  mass  removed 
and  deposited  is  described  by  the  random  variable  q.  Figure  (5.2)  illustrates  nonlocal 
mixing  in  an  idealized  one-dimensional  system  of  layers. 

It  is  further  assumed  that  a  nonlocal  mixing  event  causes  the  removal  and  de¬ 
position  of  a  fixed  amount  of  mass  -  as  opposed  to  a  fixed  mass  per  unit  volume. 
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With  this  assumption,  the  amount  of  mass  transferred  (for  nonlocal  mixing  only)  is 
not  a  function  of  the  density  at  the  location  of  removal.  The  exchange  function  for 
nonlocal  mixing  will  be  factorized  as 

p(r')Kni(r', r;  t)  =  qh( r'  -  a)p(r  -  b )q(t  -  £)  (5.25) 

and 

p(r)Kni(r, r'j  t )  =  qg( r'  -  b)h(r  -  a )q(t  -  ()  .  (5.26) 


The  shape  functions  are  defined  such  that 


(5.27) 


therefore,  they  have  units  of  inverse  volume. 

The  function  q  defines  the  time  dependence  of  the  source  and  sink  events.  If  the 
duration  of  each  event  is  small  and  there  are  many  events  occurring  within  the  time 
period  of  interest,  the  time  dependence  can  be  modeled  as  impulsive: 


q(t  -  C)  =  S(t  -  C)  •  (5-28) 

In  general,  £  and  (  are  not  independent.  For  example,  the  time  of  the  sink  event 
must  come  before  the  time  of  the  source  event,  and  the  spacing  between  events  is 
likely  small  compared  to  the  time  intervals  of  interest  in  here  (days  to  weeks). 

Using  (5.25),  (5.26)  in  (5.23),  and  grouping  all  the  nonlocal  terms  into  a  single 
forcing  term  on  the  right  hand  side  of  the  equation,  the  rate  of  change  of  mass  at 
position  r  in  the  sediment  is 

J^>(r,*)~y*  [p(r',t)Ki(r',r;t)  +  p(r,t)K,(r,r';t)]dr'  =  f(r,t)  ,  (5.29) 

where 

f(r,  i)  =  qf  [h( r'  -  a)9(r  -  b)<S(f  -  f )  -  A(r  -  a)ff(r'  -  b)<5((  -  0]  *'  .  (5.30) 
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The  left  hand  side  of  (5.29)  has  been  shown  in  Section  5.2.1  to  reduce  to  the  biodif¬ 
fusion  equation  if  K\  is  symmetric  and  local.  Two-scale  mixing  can  then  be  modeled 
as  a  biodiffusional  process  where  nonlocal  mixing  acts  as  a  forcing  function: 

A 

Wtp(T,t)  -  V  ■  [D»(r,<)V^(r,i)]  =  /(r,  t)  .  (5.31) 

Using  (5.27),  the  nonlocal  forcing  function  (5.30)  is  rewritten  as 

/( r,  i)  =  9[a(r  -  b)i((  -  ()  -  h( r  -  a )*(t  -  0]  .  (5.32) 

This  form  of  the  diffusive  mixing  equation  is  inhomogeneous  in  />;  that  is,  the  forcing 
function  is  not  dependent  on  the  value  of  p. 

Generalizing  (5.25)  and  (5.26)  as  a  series  of  source/sink  pair  functions  occurring 
randomly  in  space  and  time,  the  exchange  function  becomes 

p( r';  t)Kni(r',  r;  t)  =  ^  qnh„( r'  -  an)gn(r  -  b„)£(*  -  £„)  (5.33) 

n 

and 

p(r;  t)Kni( r,  r';  t)  =  ^  r  -  an)gn(r '  -  bn)^(<  -  Cn)  •  (5.34) 

n 

The  nonlocal  forcing  function  is  then  described  as 

/( r,  t)  =  Y^Qn  [gn{r  -  b„)5(t  -  £„)  -  hn( r  -  a n)8(t  -  Cn)]  •  (5.35) 

n 

The  source  and  sink  functions  gn  and  hn  are  random  functions  associated  with  the 
nth  event.  They  describe  a  random  shape  of  the  sediment  volume  moved  nonlocally 
by  each  event. 

In  general,  the  excitation  may  be  anisotropic  and  a  function  of  position  in  the 
sediment  volume.  For  example,  the  shape  of  an  event  may  be  asymmetric  in  the 
horizontal  and  vertical  directions,  and  biological  activity  may  decrease  with  depth 
into  the  sediment.  Figure  5.3  shows  examples  of  shape  function  that  represent:  (a) 
head-down  deposit  feeding;  (b)  burrow  hole  infilling;  and  (c)  arbitrary  movements  of 
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Figure  5.3:  Examples  of  cylindrical  and  spheroidal  source/sink  shape  functions  that 
represent  idealized  nonlocal  mixing  of  macrofauna:  (a)  head-down  deposit  feeding; 
(b)  burrow  hole  infilling;  and  (c)  arbitrary  movements  of  sediment. 


sediment.  Figure  5.4  shows  a  hypothetical  depth-dependent  distribution  of  spherical 
shape  functions.  Activity  is  concentrated  near  the  sediment-water  interface  and 
decreases  with  depth  into  the  sediment,  but  it  is  uniform  in  the  horizontal  plane. 
The  source  and  sink  functions  may  also  be  location  dependent.  Head-down  deposit 
feeders  will  deposit  sediment  only  on  the  surface  of  the  sediment,  for  example. 

The  solution  to  the  inhomogeneous  biodiffusion  equation  can  be  expressed  as  a 
convolution  integral  using  the  appropriate  Green’s  function  for  the  diffusive  system: 

p(r,*)=  f  f  gd(r,r',t-  r)/(r',  r)drdr'  .  (5.36) 

Jo  Jv 

The  Green’s  function  (gd)  is  the  solution  to  the  diffusion  equation  with  an  impulse 
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Figure  5.4:  Spherical  shape  functions  that  are  homogeneously  distributed  in  the 
horizontal  (a)  and  nonhomogeneously  distributed  in  depth  (b).  The  Poisson  point 
density  in  space  and  time  decreases  with  depth,  concentrating  activity  near  the 
sediment- water  interface. 


excitation: 

8 


j^9d(r,  t;  r0,  t0)  -  V  •  [D(r,  t)Vgd( r,  t ;  r0,  to)]  =  S( r  -  r0)S(t  -  t0)  .  (5.37) 


The  prescribed  initial  and  boundary  conditions  are 

p(x,  y,z,t  =  t0)  =  po(x,  y,  z)  , 


dz 


p(x,y,z  =  0  ,t)  =  0  . 


(5.38) 


Assuming  D  is  independent  of  r  and  i,  the  half-space  Green’s  function  for  diffusion 

is  found,  by  introducing  an  image  source  across  the  z  axis  [42],  to  be 

c-H*/[4  (t-t0)D)  +  e- R?/[4{t-t0)D] 


gd(r,t-,r0,t0 )  = 


[4ir(t-to)D]^2 


(5.39) 


where 


(5.40) 


R2  =  (x-  x0 )2 .+  (y  -  yof  +  (z-  z0f 
Ri  =  (#  —  Xof  +  (y  —  yo)2  +  (z  +  zo)2 . 

If  half-space  effects  are  neglected  (as  will  be  done  in  later  sections),  the  image  source 
is  removed  from  the  solution,  and 

e-R2/[4(t-t0)D] 


gd(r,t;r0,t0)  = 


[47r(£  —  t0)D]3/2  ‘ 


(5.41) 
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5.3.1  Perturbation  Method 

Equations  (5.31)  and  (5.35)  describe  nonlocal  mixing  in  terms  of  an  excitation  func¬ 
tion  that  is  independent  of  the  sediment  density.  This  formulation  follows  directly 
from  the  assumption  that  an  organism  moves  a  fixed  amount  of  mass  during  non¬ 
local  transport,  independent  of  the  density  at  the  location  of  the  event.  The  same 
result  can  be  derived  using  time-dependent  perturbation  theory. 

Start  by  redefining  the  nonlocal  exchange  functions  as 

Kni(v',  r;  t)  =  ^2  un^n(r/  —  a„)gfn(r  —  bn)s(t  —  £n)  (5.42) 

n 

and 


Kni( r,  r';  t)  =  ^  vngn( r'  -  bn)h(r  -  an)s(t  -  (n)  ,  (5.43) 

n 

where  vn  are  the  exchange  function  magnitudes  that  specify  the  fraction  of  the  total 
mass  present  in  a  location  that  is  removed  or  deposited.  Mass  rate  of  change  is  still 
described  by  (5.23),  but  the  excitation  has  the  form 


/( r,  t)  =  ^2  [u„(t)#„(r  -  bn)5(t  -  £„)  -  wn(r,  t)hn( r  -  a n)5(t  -  Cn)]  ,  (5.44) 


where  the  shape  function  magnitudes  are  now  defined  as 


un(f>j  —  /  vnp{v  a.)dr  , 

Jv 

wn(r,t)  =  vnp(r,t)  /  gn(v'  -b)dr'  . 
Jv 


(5.45) 


The  forced  mixing  equation  is  no  longer  inhomogeneous  in  p\  the  forcing  function  is 
dependent  on  the  value  of  p. 

Now,  assume  that  the  fraction  of  sediment  mass  per  unit  volume  moved  by  each 
source/sink  event  is  small  compared  to  initial  density  at  that  location  such  that 


\hnVd\  <  1  and  \gnVd\  <  1  • 


(5.46) 
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The  density  can  then  be  represented  in  a  perturbation  series  as 


p  —  Po  +  Pi  +  •  •  *  5  (5-47) 

where  po  is  the  solution  of  the  biodiffusion  equation  without  nonlocal  mixing.  An 
equation  for  each  order  of  density  correction  can  be  found  by  substituting  (5.47) 
into  (5.31)  with  the  new  definition  of  the  pulse  amplitudes  and  grouping  the  terms 
by  orders  of  magnitude.  The  first-order  biodiffusion  equation  is  then  found  to  be  an 
inhomogeneous  diffusion  process: 

J^>i(r,f)  -  V  •  [D(r,t)Vpx(r,f)]  =  f(r,t )  .  (5.48) 

The  excitation  function  is  a  series  of  random  pulse  functions  with  random  ampli¬ 
tudes,  as  in  (5.44).  However,  the  event  amplitude  functions  are  now  defined  as 

^n(^)  =  /  Vnpo(r ,  t)hn(r  a„)dr  , 

Jv  f  (5.49) 

wn(r,t)  =  vnp0(r,t)  /  gn( r'  -  bn)dr'  , 

Jv 

which  are  independent  of  the  first-order  fluctuation  in  density.  Using  (5.49)  and 
comparing  (5.44)  with  (5.35),  the  fixed  mass  assumption  can  be  interpreted  in  terms 
of  the  perturbation  result.  In  doing  so,  the  mass  moved  by  each  event  is  defined  as 

qn  =  vnp0  .  (5.50) 

Either  approach  can  be  used  to  derive  the  inhomogeneous  biodiffusion  equation. 

5.3.2  Poisson  Pulse  Process 

Consider  the  excitation  of  the  inhomogeneous  biodiffusion  process.  If  the  correlation 
of  the  excitation  function  can  be  found,  the  correlation  and  corresponding  spectrum 
of  the  sediment  density  fluctuations  can  be  found.  By  assuming  the  sequence  of  exci¬ 
tation  pulse  functions  is  a  Poisson  process  in  time,  the  problem  is  greatly  simplified. 


122 


To  cast  the  problem  as  a  Poisson  process,  two  assumptions  are  made  about  the 
stochastic  properties  of  the  random  time  variables:  1)  the  random  variables  £n  and 
Cn  are  statistically  independent  and  identically  distributed;  and  2)  the  series  con¬ 
sists  of  N  pulses  that  are  homogeneously  distributed  over  a  time  interval  of  length 
T  with  N  — ►  oo  and  T  -*  oo.  The  probability  density  functions  (pdfs)  of  all  the 
random  time  variables  (£i,  •  •  •  ,  £j\r,Ci,  •  •  •  ,  Cn)  are  then  described  as  the  product  of 
the  individual  pdfs: 

p(6»Ci>”*  >€n,(n)  =  ■  (5-5i> 

n 

The  independence  assumption  essentially  amounts  to  randomizing  the  source/sink 
pair  functions  in  time.  The  sink  that  is  related  to  a  particular  source  can  occur 
at  any  time  before  or  after  the  source  event.  This  violates  causality  in  the  mixing 
process.  However,  if  a  large  number  of  events  occur  between  observation  times  and 
initial  condition  effects  are  ignored,  the  cumulative  effects  of  the  events  will  be  the 
same  as  if  their  order  was  preserved. 

Similar  assumptions  are  also  made  about  the  random  position  variables.  Up 
to  this  point,  the  statistical  properties  of  the  random  variables  an  and  b„  have 
been  general.  Now,  it  is  assumed  that  each  event  is  independent  and  identically 
distribution,  such  that 

p(ai,bi,---  ,ajy,bw)  =  JJpa&(an,bn)  .  (5.52) 

n 

L 

The  random  positions  an  and  bn  for  a  single  event  n,  however,  are  not  independent. 
In  general,  the  position  an  may  be  conditional  on  the  position  bn  for  the  same  n. 
This  may  be  the  case  for  specific  types  of  nonlocal  deposit  feeding,  such  as  head  down 
feeding  where  the  position  of  the  deposition  is  most  likely  above  the  position  of  the 
removal  of  sediment,  for  example.  It  is  further  assumed  that  the  random  position 
variables  are  independent  of  any  other  random  variables  that  may  be  defined  to 
describe  the  shape  functions  g  and  h. 
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The  correlation  of  the  random  process  /  in  space  and  time  is  defined  as 

Rff(ri,r2,tut2)  =  {f(ruU)f(r2,t2))  .  (5.53) 

Using  (5.35)  and  assuming  stationarity  in  time  only  (with  r  =  ti~t2),  the  correlation 

of  the  excitation  can  be  expressed  as  a  double  summation  with  N2  terms: 

N  N 

Rff(ri,r2,T)=^2^2[(^9nqmgmS(t-^n)S(t-T-(m)) 
m  n 

(5.54) 

Cn)&(t  T  Cm)) 

"f"  {9nhn9m^m^(t  Cn)&(t  T  Cm))]  • 

The  angle  brackets  are  the  ensemble  averages  over  the  random  variables  an,  bn,  Cn, 
Cn,  and  any  other  random  parameters  that  define  the  shape  functions  (yet  to  be 
defined).  The  arguments  to  the  shape  functions  gn  and  hn  have  been  dropped  to 
save  space. 

Using  the  Poisson  pulse  assumption,  the  random  time  variables  with  m^n  are 
assumed  to  be  independent.  The  double  summation  can  then  be  divided  into  two 
parts: 

N 

Rffi rl,  r2,  r)  =  [{qngnqn9n){S(t  -  £n)6(t  -  T  -  Cn)) 

n 

2(^n^nQn/ln)(5(t  T  Cn)) 

+  ( qnhnqnhn)(S(t  -  Cn)S{t  T  Cn))] 


N 


(5.55) 


+  Y  l(<h9n)(qm9m)(S(t  ~  £»))($(*  -  T  -  Cm)) 

m^n 


-  2{qngn){qmhm)(5(t  -  Cn))(S(t  -  r  -  Cm)) 

+  (qnhn)(qmhm){S(t  -  Cn)){S(t  -  T  -  Cm))]  ■ 

The  function  ( qngn )  is  the  first  moment  of  the  source  function  (to  be  discussed  later), 
defined  as 


{qngn)  =  {qngn(  r  -  bn)) 


(5.56) 


/ 
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The  function  ( qnhn )  is  the  first  moment  of  the  sink.  The  function  (qngnqngn)  is  the 
second  moment  of  the  source  function,  defined  as 


( qhQnqnQn )  —  bji)*7n( ^2  hn))  • 


(5.57) 


Similarly  for  ( qngnqnhn )  and  {qnhnqnhn). 

Since  the  random  time  functions  are  impulsive  and  homogeneously  distributed 
over  time  T,  the  ensemble  averaging  over  time  can  be  performed,  resulting  in 


m  -  C»W*  -  T  -  in))  =  fS[r)  ■ 


(5.58) 


In  addition,  since  the  random  variables  are  independent  and  identically  distributed, 
the  summations  can  be  eliminated.  The  first  summation  contains  N  identical  terms, 
and  the  second  summation  contains  N(N  —  1)  identical  terms.  Equation  (5.55)  is 
then  rewritten  as 


Rf  }{r  i,r  2,  r) 


at  9N 

=  2t[(92^)  +  (g2hh)]5(T)  -  j^(q2gh) 

+  n{nt7  [(«>  -  w]2  • 


(5.59) 


The  angle  brackets  now  represent  the  ensemble  average  over  the  random  variables 
a,  b  and  the  other  random  parameters  that  define  g  and  h.  As  N  and  T  become 
large,  the  Possion  point  density  in  time  is  defined  as 


(5.60) 


The  second  term  in  (5.59)  will  then  vanish  in  the  limit,  and  the  excitation  correlation 
can  be  rewritten  as 


i?//(ri,r2,r)  =  Xt[(q2gg)  +  ( q2hh)]8(r )  +  X2[(qg)  -  (< qh )]2  (5.61) 


In  general,  Xt  can  be  a  function  of  time.  This  is  known  as  an  inhomogeneous  Poisson 
point  process  in  time. 
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To  solve  for  the  correlation  of  the  excitation  function,  the  correlations  of  the 
shape  functions  must  be  found.  Before  any  assumptions  are  made,  it  is  useful 
to  describe  the  random  shape  functions  as  deterministic  three-dimensional  pulse 
functions  that  are  functions  of  a  finite  number  of  random  parameters: 


h( r  —  a)  =  h( c,  r  —  a)  , 
flf(r-b)  =  $(d,r-a)  . 


(5.62) 


The  randomness  in  the  source  and  sink  functions  ( g  and  h)  is  now  represented  by 
the  random  variables  c  and  d.  In  general,  c  and  d  can  be  random  vectors  of  a 

A 

finite  number  of  random  variables  each.  The  deterministic  functions  g  and  h  now 
represent  the  nonrandom  characteristic  shape  of  the  source  and  sink  events. 

In  general,  the  magnitudes  of  the  shape  functions  are  not  independent  of  the 
parameters  that  define  the  randomness  in  the  shape  function.  For  example,  the 
variable  c  may  represent  the  random  radius  of  the  sink  function  h.  The  magnitude 
of  the  sink  event  ( q )  could  also  be  a  function  of  the  random  radius  -  a  larger  volume 
of  sediment  affected  will  result  in  a  larger  amount  of  mass  that  is  transported.  If 
this  is  the  case,  the  joint  probability  density  function  for  q  and  c  must  be  specified 
to  find  the  correlation  of  the  sink  event. 

However,  if  the  magnitudes  of  the  source/sink  events  are  assumed  to  be  inde¬ 
pendent  of  the  shape  functions,  the  correlations  can  be  simplified  as 


(gg)  =  (9>(<Kd>r-b))  , 

(g2gg)  =  (g2)(g{ -  b)£(d,r2  -  b)) . 


(5.63) 


The  correlation  function  can  then  be  described  in  terms  of  the  correlations  of  the 
shape  functions  only: 


#//(r1,r2,r)  =  V92[i?as(n,r2)  +  iMri,r2)]£(r)  +  *Wq[{g)  ~  ( h )f  »  (5-64) 

where  cr2  is  the  variance  of  the  event  magnitudes  and  rjq  is  the  mean.  The  shape 
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function  correlations  are  defined  as 

Rgg{ ri,r2)  =  {g(d,r!  -  b)^(d,r2  -  b))  , 

Rhh{ ri,r2)  =  (fi(c,ri  -  a)/>(c,r2  -  a))  . 

The  covariance  of  the  excitation  function  is  defined  as 

C}f{rur2)  =  Rff( ri,r2)  -  ( f(r,t ))2  . 

It  can  be  easily  shown  that  the  second  term  in  (5.64)  is  the  equal  to  the  second  term 
in  (5.66).  Therefore,  the  covariance  of  the  excitation  is 

C//(ri,r2,r)  =  A^[f?aff(rx,r2)  +  Rhh{ru r2)]<S(r)  .  (5.67) 

5.3.3  Stationarity  in  Space 

Up  to  this  point,  the  only  assumptions  made  about  the  statistics  of  the  shape  func¬ 
tions  in  space  are  that  the  source/sink  events  are  independent  and  identically  dis¬ 
tributed.  The  pdfs  of  the  random  variables  that  define  the  shape  functions  are  still 
general.  The  correlations  of  the  shape  functions  could  have  arbitrary  structure  that 
models  the  nonstationary  behavior  of  complicated  systems  of  organisms.  Variability 
in  the  horizontal  dimension  due  to  the  clustering  of  organisms  could  be  modeled  as  a 
random  process  with  stationary  increments,  for  example.  The  depth  dependence  of 
bioturbation  could  be  modeled  by  distributing  the  depth  of  the  source/sink  events 
nonuniformly,  concentrating  activity  near  the  sediment-water  interface.  The  pdf  of 
the  random  position  variables  a  and  b  may  also  be  joint  and/or  conditional  upon 
each  other.  In  general,  some  form  of  non-stationarity  in  the  depth  direction  is  ex¬ 
pected,  unless  the  distribution  of  source/sink  events  is  uniform  in  depth.  Uniformity 
in  depth  is  unlikely  and  is  a  very  restrictive  assumption. 

This  inherent  non-stationarity  presents  an  apparent  paradox  in  relating  a  model 
for  bioturbation  to  the  acoustic  model  of  Chapter  2.  Somehow,  the  spectrum  of 
the  bioturbational  process  must  be  found,  and  a  spectrum  is  only  defined  in  terms 


(5.65) 


(5.66) 
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of  a  stationary  process.  Alternatively,  the  acoustic  correlation  could  be  found  by 
integrating  the  sediment  correlation  directly.  However,  this  complicates  the  inverse 
process  of  inferring  biological  parameters  and  forfeits  the  physical  insight  that  can 
be  gained  by  working  with  spectra. 

This  apparent  paradox  in  modeling  can  be  resolved  by  introducing  a  locally 
stationary  process.  This  leads  to  the  concept  of  a  local  spectrum  that  can  vary  in 
space.  The  simplest  form  of  a  local  spectrum  is  found  by  assuming  the  positions 
of  the  source/sink  events  are  independent  and  uniformly  distributed  with  a  Poisson 
density  that  is  a  function  of  position  in  space.  The  shape  function  correlations  can 
then  be  defined  in  terms  of  the  difference  coordinate  =  ri  —  r2,  such  that 


Rhh{rd)  =  JJ  A0(a)h(c,  r  -  a)A(c,  r  -  rd  -  a)pc(c)  dc  da 
=  JJ  A0(r  -  a)h(c,a)h(c,a  —  rd)pc(c)  dc  da  , 


(5.68) 


where  pc(c)  is  the  pdf  of  the  random  variable  c,  r  =  rx,  and  the  change  of  variable 
a  =  r  —  a  is  made.  The  parameter  Aa  is  the  nonhomogeneous  Poisson  density  in 
space.  If  A0  varies  slowly  over  the  integral  on  a ,  compared  to  the  shape  function  h, 
the  Poisson  parameter  can  be  moved  outside  of  the  integral: 

Rhh(rd)  =  Aa(r)  //Me,  a)h( c,  a  —  rd)pc(c)  dc  da  .  (5.69) 

Similarly  for  the  correlation  function  Rgg : 

R99(rd)  =  Afc(r)  JJ g(d,a)g(d,a  -  rd)pd(d)  dd  da  .  (5.70) 

The  excitation  correlation  is  now  simplified  as  a  position-dependent,  locally  station¬ 


ary  process, 


where 


C/f(rd,r,T)  =  a2q[\b{r)Rgg{rd)  +  Xa(r)Rhh(rd)]5(T)  , 


(5.71) 


Rgg(rd )  =  JJ g(d,a)g(d,a-  rd)pd(d)  dd  da  (5.72) 
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and 


Rhh{rd)  =  ilk  c,a)h(c,a  —  rv)pc(c)  dc  da  .  (5.73) 

The  Poisson  parameter  in  time  has  been  absorbed  into  Xa  and  A&,  which  are  redefined 
as  the  number  of  source  events  (Nb)  or  the  number  of  sink  events  (Na  )  per  unit  time 
and  unit  volume: 


\  ^0.  ,  y.  Nb 

Xa  =  —  and  Xb  = 


TV 


TV 


(5.74) 


5.4  Spectrum  of  a  Biodiffusional  Process 


Equation  (5.31)  is  a  deterministic  differential  equation  with  a  random  excitation  / 
due  to  the  nonlocal  exchange  of  sediment.  The  solution  p  can  be  interpreted  as  the 
output  of  a  linear  system  with  a  stochastic  input: 


p(r,t)  =  Ld[f(r,t)}  .  (5.75) 

The  operator  Ld  is  the  integral  equation  (5.36)  containing  the  diffusion  Green’s 
function.  If  the  half-space  effects  are  neglected,  and  the  diffusion  process  is  linear 
and  space-  and  time-invariant,  p  can  be  expressed  as  a  space-time  convolution  with 
the  system  impulse  response, 

p(r,  t )  =  gd{ r,  t )  ®  /( r,  t)  ,  (5.76) 


where 

a,(r,<)  =  Ij[f(rM(l)]  •  (5.77) 

It  follows  from  the  linearity  of  the  expected  value  that  the  correlation  of  the  output 
is  found  in  terms  of  the  correlation  of  the  input  process: 

RPP(rd,  t)  =  Rf/(rd,  r)  ®  gd(rd,  r)  ®  gd(- rd,  -r)  .  (5.78) 
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Thus,  the  mixing  process  can  be  characterized  as  a  linear  filtering  process  where 
the  spectrum  of  the  output  process  ( Spp )  is  defined  completely  in  terms  of  the  spec¬ 
trum  of  the  input  (S/j)  and  the  frequency  response  of  the  diffusive  system  (Gd): 

S„( k,«)  =  |GJ(k,w)|2S//(k,u,)  .  (5.79) 

The  function  Gd  is  the  space-time  Fourier  transform  of  gd .  The  spectrum  of  the 
input  is  found  by  taking  the  Fourier  transform  of  the  locally  stationary  covariance 

function  (5.71)  with  respect  to  r<f  and  r, 

2 

.S//(k, r)  =  g-  [A»(r)S„(k)  +  A„(r)SM(k)]  ,  (5.80) 

where 

Shh(  k)  =  (27r)3(|7f(k)|2)  (5.81) 


and 


S93(  k)  =  (27r)3<|G(k)|2)  (5.82) 


are  the  spectra  of  the  shape  functions,  and  G  and  H  are  the  Fourier  transforms  of 
the  shape  functions,  g  and  h.  The  angle  brackets  in  (5.81)  and  (5.82)  are  now  the 
expected  value  operators  with  respect  to  the  random  variables  c  and  d  only.  Note 
that  the  spectrum  is  independent  of  temporal  frequency,  w,  because  of  the  assumed 
impulsive  nature  of  the  source  and  sink  events. 

For  the  case  of  a  depth-independent  and  isotropic  diffusion  coefficient  and  ne¬ 
glecting  the  half-space,  the  diffusion  Green’s  function  is 

<«*•">  "Z <5-83> 


where  k2  =  +  ky  +  k*.  The  spectrum  Spp  is  then  expressed  as 


Spp(k,  r,w)  = 


SffjK  r) 
T>2fc4  +  u;2  ‘ 


(5.84) 
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Performing  an  inverse  Fourier  transform  in  time,  the  temporal  and  spatial  spectrum 
of  density  fluctuations  in  the  sediment  can  be  modeled  as 

S„(  k,r,r)=^^e-D‘’M.  (5.85) 

5.5  Inhomogeneous  Biodiffusion  of  Sediment  Microtopography 

Epifaunal  and  infaunal  activity  create  time  varying  roughness  at  the  sediment-water 
interface  on  spatial  and  temporal  scales  that  are  of  interest  in  acoustic  scattering. 
Mega-  and  macrofaunal  locomotion  along  the  interface  create  surficial  traces.  To¬ 
pographical  trace  concentration  has  been  considered  as  an  indicator  of  abundance 
of  epifauna  [84].  Infauna,  such  as  deposit  feeders,  also  create  and  destroy  microto¬ 
pography  as  sediment  is  transported  to  and  from  the  sediment  interface,  creating 
mounds  and  depressions  on  the  surface. 

While  larger  organisms  create  roughness,  other  processes  act  to  erase  their  effects. 
Physical  processes,  such  as  gravity  and  fluid  flow  due  to  wave  action  or  currents,  will 
transport  surface  sediment  particles  and  erase  traces  of  larger  animals.  Meifaunal 
activity  (analogous  to  the  activity  on  similar  scales  within  the  sediment)  creates 
random  and  continuous  local  exchange  at  the  surface  acting  to  erase  or  smooth 
larger-scale  roughness  features. 

Analogous  to  the  two-scale  mixing  of  sediment  volume  inhomogeneities,  as  de¬ 
tailed  in  Section  5.3,  the  temporal  and  spatial  properties  of  surface  roughness,  due  to 
reworking  by  benthic  organisms,  can  also  be  modeled  as  a  forced  diffusion  process. 
The  methods  used  to  describe  the  correlation  and  spectral  properties  in  terms  of 
animal  activity  and  spatial  characteristics  are  applied  directly  to  the  surface  rough¬ 
ness.  The  diffusive  quantity  is  now  the  rough  surface  height  (zs),  and  local  mixing  by 
meiofauna  (and  other  physical  mechanisms)  is  described  with  a  horizontal  diffusion 
coefficient  for  roughness  ( Ds ). 

Equation  (5.31)  for  two-scale  mixing  of  sediment  density  inhomogeneities  is  then 
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rewritten  to  describe  the  temporal  and  spatial  evolution  of  surface  relief: 

JU.(r,f)  -  V  ■  [B»(r,t)Vz,(r,«)]  =  /.(r,t)  .  (5.86) 

The  position  vectors  and  derivatives  are  now  in  two  dimensions  r  =  (x,y).  The 
nonlocal  forcing  function  is  described  as 

fs(r,t)  =  Y^<lnhn(r-an)S{t-in)  ■  (5.87) 

n 

The  magnitude  of  each  event  is  defined  by  the  zero-mean  quantity  qn  in  units  of 
height.  The  shape  function  hn  is  a  dimensionless  random  function  associated  with 
the  nth  event  that  describes  a  random  shape  of  a  trace  left  behind  by  an  animal 
or  other  surface  relief  event.  The  diffusion  coefficient  Ds  describes  the  localized 
mixing  in  the  horizontal  direction  at  the  surface  of  the  sediment.  In  some  cases,  it 
may  be  the  same  as  the  horizontal  component  of  the  biodiffusion  coefficient  Db  used 
in  Section  5.3.  However,  it  may  also  may  include  processes  that  only  exist  on  the 
sediment-water  interface,  such  as  mixing  due  to  bottom  currents. 

The  correlation  of  the  excitation  function  is  found  in  the  same  manner  as  it  was 
found  in  Section  5.3.2.  However,  now  it  is  the  the  covariance  of  fs  that  is  defined 
in  terms  of  the  correlation  of  the  shape  functions,  the  variance  of  the  magnitudes  of 
the  events,  and  the  the  Poisson  density: 

Cff( ri, r2, r)  =  r2)5(r)  .  (5.88) 

The  variance  of  the  event  magnitudes  qn  is  cr2.  The  Poisson  density  As  is  the  number 
of  surficial  traces  per  unit  time. 

The  problems  associated  with  stationarity  in  the  depth  direction,  as  in  the  volume 
problem,  do  not  exist  in  the  two-dimensional  surface  roughness  problem.  Stationar¬ 
ity  in  the  horizontal  directions  is  more  easily  justified  (on  spatial  scales  comparable 
with  the  sizes  of  the  traces  left  by  mega-  and  macrofauna,  at  least).  Therefore, 
with  \s  defined  as  a  slowly  varying  function  of  space  compared  with  the  trace  shape 
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functions,  the  covariance  of  the  excitation  function  can  be  modeled  using  Section 
5.3.3,  and  the  spectrum  of  the  excitation  modeled  using  Section  5.4,  such  that 

S//(k,r)  =  |U.(r)S»*(k) .  (5.89) 

The  Poisson  density  As  is  now  the  number  of  animal  trace  events  per  unit  time 
per  unit  surface  area,  and  Shh  is  the  spectrum  of  the  trace  shape  functions  defined 
analogously  to  (5.81).  The  surface  roughness  spectrum  of  a  biodiffusive  process  at 
the  sediment-water  interface  is  then 


Szz(k,r,a>) 


%(k>r) 


(5.90) 


The  temporal  and  spatial  spectrum  of  surface  roughness  due  to  biological  mixing  is 


&*(k,  r,  r) 


S//(k» r)  -P,fc2l rl 

2  Dsfc2 


(5.91) 


5.6  Summary  and  Conclusions 


The  above  expressions  are  general  in  the  sense  that  they  attempt  to  model  biotur- 
bation  in  nonspecific  terms  (e.g.,  diffusion  and  point  density)  while  still  accounting 
for  the  two-scale  nature  of  the  mixing  process  (meiofauna  and  macrofauna)  and 
non-stationarity.  The  shape  functions  and  their  expected  values  should  reflect  the 
size  and  shape  of  the  macrofauna  that  move  sediment  nonlocally.  Anisotropy  can  be 
modeled,  for  example,  by  using  a  shape  function  that  represents  directional  feeding, 
such  as  vertical  burrows  made  by  conveyer-belt  feeders  [40,  7].  Spatial  dependence 
of  the  source/sink  events  is  modeled  as  a  nonhomogeneous  Poisson  process  in  space 
and  time.  The  preference  of  a  deposit  feeder  to  remove  sediment  from  a  specific 
depth,  for  example,  could  be  modeled  by  specifying  a  non-uniform  Poisson  density 
function  that  increases  the  probability  of  an  event  occurring  at  the  favored  depth. 

The  Poisson  point  density  in  time  should  reflect  the  amount  and  rate  of  macro- 
faunal  activity,  which  in  general  can  be  time  dependent.  The  diffusion  coefficient 
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should  reflect  the  length  and  time  scales  of  the  meiofaunal  activity,  which  can  also 
be  modeled  as  depth  dependent.  Of  course,  in  some  circumstances  the  Poisson  as¬ 
sumption  (homogeneous  or  nonhomogeneous)  cannot  be  justified.  This  may  occur 
if  the  gradient  in  the  Poisson  parameter  cannot  be  assumed  constant  over  the  ge¬ 
ometric  extent  of  the  shape  function.  For  example,  deposit  feeders  with  long  body 
lengths  can  extend  deep  into  the  sediment.  In  this  case,  equation  (5.67)  must  be 
solved  directly,  as  opposed  to  using  (5.71).  The  possibility  of  describing  mixing  as 
a  stationarity  process  may  be  abandoned  altogether.  This  would  bring  in  question 
the  use  of  spectral  analysis  for  describing  sediment  physical  properties.  Correlation 
descriptions  would  still  be  valid,  however,  using  (5.67). 
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Chapter  6 

REMOTE  SENSING  OF  BENTHIC  ACTIVITY 

Remote  sensing  of  benthic  biology  by  acoustic  wave  scattering  is  potentially  an 
important  tool  for  studying  biological  processes  on  spatial  and  temporal  scales  larger 
than  can  be  achieved  by  current  point  sampling  methods.  Inferences  about  the  ef¬ 
fect  of  benthic  biota  on  sediment  structure  are  typically  made  by  examining  sparsely 
sampled  sediment  cores  and  optical  imagery.  Coring  can  provide  detailed  informa¬ 
tion  on  the  vertical  structure  of  the  sediment  at  an  instant  in  time  and  for  a  single 
horizontal  location.  However,  large  gaps  are  generally  left  in  horizontal  and  temporal 
sampling.  High  resolution  X-radiographs  [28],  X-ray  tomography  [52],  microelectric 
conductivity  methods  [74],  and  box  coring  provide  better  horizontal  resolution,  but 
they  are  usually  limited  to  very  small  spatial  scales  (less  than  a  meter).  They 
are  also  destructive,  providing  only  a  single  measurement  in  time.  Nondestructive 
in  situ  methods,  such  as  acoustic  tomography  in  the  sediment  [74,  87]  and  very 
high-frequency  acoustic  (ultra-sound)  imaging  [51],  can  provide  increased  temporal 
resolution  with  minimal  impact  on  the  sediment.  However,  they  are  also  limited 
to  small  spatial  scales  due  to  high  acoustic  attenuation  in  the  medium.  Optical 
methods,  such  as  stereo  photography  [64,  70],  only  provide  information  about  the 
sediment  interface,  and  are  again  limited  in  their  larger-scale  spatial  resolution. 

This  chapter  will  discuss  a  stochastic  model-based  technique  for  monitoring  bi¬ 
ological  activity  on  large  spatial  and  temporal  scales  using  high-frequency  acoustic 
backscatter.  The  method  is  stochastic  because  it  relies  on  the  random  nature  of 
biological  mixing  to  create  the  scattering  medium.  In  this  sense  it  is  similar  to  pre¬ 
vious  statistical  methods  of  monitoring  benthic  activity  using  backscatter  [41].  The 
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method  differs  in  its  reliance  on  physical  models  of  scattering  and  sediment  mixing. 

The  chapter  will  begin  with  the  development  of  a  forward  model  of  scattering  due 
to  bioturbation  using  the  results  of  chapters  2  and  5.  Then  the  inverse  problem  will 
be  investigated  where  an  estimate  of  the  temporal  correlation  of  backscatter  is  used 
to  infer  parameters  of  biological  mixing  in  the  sediment.  Finally,  the  models  will 
be  applied  to  experimental  data  as  a  preliminary  test  of  the  modeling  methodology. 
The  experiment  was  performed  in  the  West  Sound  of  Orcas  Island  in  Puget  Sound. 
Backscatter  data  was  recorded  in  an  area  with  observed  biological  activity.  Ground 
truth  estimates  of  biological  activity  are  compared  with  the  acoustically  inferred 
activity. 

6.1  Forward  Model  of  Scattering 

A  model  for  the  effects  of  bioturbation  on  volume  scattering  is  developed  by  relating 
the  decorrelation  in  time  of  the  backscattered  field  with  changes  in  the  sediment 
physical  properties  due  to  biological  mixing.  The  perturbation  method  of  Section 
2.4  is  used  to  describe  scattering  as  a  function  of  the  statistical  properties  of  the 
sediment.  The  sediment  bioturbation  model  of  Chapter  5  is  used  as  the  spectral 
model  of  sediment  mixing.  They  are  combined  to  form  a  forward  stochastic  model 
of  scattering  due  to  benthic  biological  activity. 

Consider  backscatter  of  an  acoustic  field  due  to  volume  scattering  in  sediment. 
For  backscatter  at  grazing  angle  0,  equations  (2.74)  and  (2.56)  are  rewritten  as 

cpp(e,T)=s^\mi2mt.+if  (6.i) 

where 

iW  =  -7r%e‘lh,r .  (6.2) 

p  (47T r)z 

and  kd  =  (Re[&]/2)  sin  6  is  the  Bragg  wavenumber.  Recall  that  Spp(kd,T )  is  the 
time-dependent,  half-space  spectrum  of  density  fluctuation  in  the  sediment  defined 
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using  (2.75)  as 

5„(k*r)  =  f  drc  f°°  dk'dz  SJkdx,kdy,k'dz,rc,T)W(k'dz  -  kdz,zc)e~^  ,  (6.3) 

J  hs  J — oo 

where 


W(kz,zc) 


sin  (2  kzzc) 
nkz 


(6.4) 


The  function  5Pp(krf,rc,r)  is  the  time  and  space  dependent  spectrum  of  density 
fluctuations  in  the  sediment  evaluated  at  the  Bragg  wavenumber.  The  normalized 
correlation  coefficient  for  scattering  is  defined  as 

Gpp{0 ?  T )  $ ppifedi  T*) 


Cpp(T )  “ 


(6.5) 


Cpp(0, 0)  Spp(kd,  0) 

In  Section  5.4  a  model  for  the  spectrum  of  density  fluctuations  due  to  biological 
mixing  was  presented,  see  equations  (5.85)  and  (5.80).  Mixing  is  modeled  as  a  locally 
stationary  and  time-dependent  diffusion  process  having  the  spectrum 


q  /l-  r  T\  _  v)  -d&\t\ 

bpp{\z,r,T)-  wk2 


(6.6) 


where  the  spectrum  of  the  random  biological  forcing  is 


Sff( k,r)  =  47TVJ  [A6(r)(|(7(k)|2)  +  Aa(r)(|i7(k)|2)]  .  (6.7) 

Equations  (6. 1-6.7)  form  the  basis  of  the  forward  and  inverse  models  of  scattering 
due  to  biological  activity.  The  model  parameters  include:  the  biodiffusion  coefficient 
D;  the  spectra  of  the  macrofaunal  shape  functions  h  and  g\  and  the  nonhomogeneous 
Poisson  point  density  functions  Aa  and  A &.  The  observable  is  the  estimate  of  the 
correlation  function  (6.1). 

Using  (6.1),  the  correlation  of  the  scattered  field  is  a  function  of  lag  time  and 
grazing  angle.  As  the  lag  time  increases,  the  correlation  is  expected  to  decrease 
exponentially  with  a  time  constant  that  is  predicted  by  the  time-dependent  biodif- 
fusive  spectrum  of  density  fluctuations  using  (6.6).  Therefore,  the  time  dependence 
of  the  field  correlation  is  due  to  the  time  dependence  of  the  diffusion  process. 
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The  grazing  angle  dependence  of  the  backscatter  correlation  is  due  to  the  scat¬ 
tering  geometry.  At  low  grazing  angles,  the  incident  field  will  penetrate  a  short 
distance  into  the  sediment.  Therefore,  the  scattered  field  at  low  grazing  angles  will 
be  influenced  primarily  by  the  sediment  near  the  sediment-water  interface.  At  higher 
grazing  angles,  the  incident  field  will  penetrate  deeper  into  the  sediment,  and  the 
scattered  field  will  be  affected  by  deeper  sediment,  as  well  as  the  interface  sediment. 
If  the  sediment  mixing  processes  are  depth  dependent,  the  grazing  angle  depen¬ 
dence  of  the  scattered  field  correlation  should  reveal  this  structure.  A  concentration 
of  biological  activity  near  the  sediment-water  interface  (decreasing  with  depth)  will 
cause  the  scattered  field  at  lower  grazing  angles  to  decorrelate  faster  than  at  higher 
grazing  angles. 

6.2  Inverse  Model  of  Scattering 

Using  the  forward  model  of  (6.1)  through  (6.7),  predictions  of  backscatter  corre¬ 
lation  can  be  made  using  information  about  the  parameters  that  define  biological 
mixing.  If  the  stochastic  geometry,  behavior,  and  rate  of  activity  of  an  organism 
are  well  understood,  its  effect  on  the  scattered  field  over  time  can  be  modeled.  In 
contrast,  the  remote  sensing  problem  is  the  estimation  of  these  stochastic  biological 
parameters  from  the  observation  of  the  decorrelation  of  acoustic  backscatter. 

6.2.1  Estimation  of  Backscatter  Correlation 

To  use  the  backscattered  field  as  a  tool  for  observing  changes  in  the  sediment,  the 
inverse  problem  must  begin  by  estimating  the  temporal  correlation  function  (6.1) 
using  acoustic  data.  One  means  of  doing  this  is  to  use  the  time-domain  scattered 
field  to  form  an  estimate  of  the  ensemble  average.  To  start,  a  finite  duration  pulse 
is  transmitted  at  time  to ,  and  the  time-domain  scattered  field  is  digitized  and  time¬ 
gated  to  correspond  to  scattering  from  an  isolated  and  finite  region  of  the  seafloor. 


138 


An  estimate  of  the  backscatter  correlation  can  be  found  by  correlating  the  time- 
domain  scattered  field  from  a  volume  of  sediment  the  size  of  the  ensonfied  volume. 
If  the  ensonified  volume  is  large  compared  to  the  correlation  length  of  the  medium, 
the  time-domain  correlation  will  be  an  estimate  of  the  correlation  function  (2.74) 
(see  Section  2.6). 

A  pulse  is  transmitted,  and  the  received  signal  ps(t)  is  digitized  as 

=  ,  (6.8) 

with 

«.(<„)  =  .  (6.9) 

The  discretized  variable  tn  is  the  time  after  transmission  of  the  pulse  for  a  single 
scan.  The  backscattered  signal  from  each  range  bin  is  then  cross-correlated  with  the 
signal  from  the  same  range  bin,  but  from  a  scan  at  some  later  time  r: 

M+N- 1 

A>p(r)=  X  «•(*»)«:(*»+ *-)  •  (6-10) 

n—M 

The  signal  us(tn)  is  from  the  scan  chosen  as  a  reference,  and  ua(tn  +  r)  is  from  a 
scan  at  time  r  later.  The  limits  of  the  summation  (M  and  N)  are  the  start  and  end 
sample  numbers  for  the  range  bin  of  interest. 

The  correlation  estimate  is  a  measure  of  the  total  change  in  the  scattered  field 
between  scans,  including  changes  in  the  propagation  path  between  the  transducer 
and  the  seafloor.  To  be  assured  that  any  decorrelation  observed  in  the  scattered 
field  is  due  to  changes  in  the  seafloor  only,  propagation  effects  must  be  removed. 
The  most  likely  source  of  decorrelation  introduced  by  the  propagation  path  is  sound 
speed  changes  in  the  water.  Uncorrelated  ambient  noise  and  volume  reverberation 
from  the  water  due  to  suspended  matter  and  organisms  [41],  for  example,  are  other 
sources  of  signal  decorrelation. 
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The  speed  of  sound  in  the  ocean  water  will  fluctuate  with  changes  in  the  wa¬ 
ter  temperature.  Neglecting  the  scattering  from  temperature  microstructure  in  the 
water,  changes  in  the  average  water  temperature  between  scans  of  the  seafloor  will 
introduce  a  phase  and  amplitude  shift  into  the  scattered  field,  as  compared  with 
the  reference  field  [20].  Before  cross-correlating  the  scattered  fields  from  different 
scans,  the  sound  speed  fluctuations  due  to  water  temperature  changes  (Sc)  must  be 
removed.  The  temperature  compensated  cross-correlation  between  scans,  derived 
previously  in  [20],  can  be  written  as 

M+N-l 

C„(t)  =  Y.  «.(*")<(*»  +  .  (6.11) 

n=M 

This  approximation  compensates  for  the  shift  in  the  carrier  (at  frequency  u>o),  but 
neglects  the  changes  in  the  complex  envelope  of  the  scattered  field. 

6.2.2  Models  of  Biological  Activity 

The  mixing  model  of  (6.6)  is  very  general.  With  no  a  priori  information  about  the 
statistical  description  of  the  organisms  in  a  sediment  to  constrain  the  spatial  and 
temporal  parameters  of  the  bioturbation  model,  the  remote  sensing  problem  may  be 
ill-posed.  It  is  likely  that  the  spectrum  of  the  shape  functions  will  provide  the  most 
ambiguity  in  the  model  inversion,  as  the  types  of  shape  functions  (that  reflect  the 
types  of  organisms)  and  their  spatial  distributions  provide  the  most  unconstrained 
degrees  of  freedom  in  the  model.  It  is  unrealistic  to  expect  a  unique  inversion  without 
some  prior  knowledge  of  their  form.  In  general,  some  knowledge  of  the  biology  must 
be  applied  to  the  inverse  problem. 

A  formal  analysis  of  the  stability  of  the  inversion  will  not  be  performed.  However, 
several  assumptions  about  the  sediment  source/sink  functions  are  made  to  simplify 
the  sediment  mixing  model,  and  to  improve  the  uniqueness  of  the  inversion.  First, 
the  biodiffusion  coefficient  is  assumed  to  be  uniform  in  the  horizontal  directions  and 
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exponentially  decaying  in  depth: 

D(z)  =  D0e~WL  .  (6.12) 

This  form  assumes  that  diffusive  activity  decreases  with  depth  at  a  rate  defined 
by  the  meiofauna  rework  depth  L ,  with  meiofaunal  activity  concentrated  near  the 
surface.  It  also  seems  reasonable  to  assume  the  same  form  for  the  spatial  point 
density  (the  Poisson  parameters  Aa  and  A&): 

A.(r)  =  A0e-|2|/L*  .  .  (6.13) 


The  activity  is  assumed  to  decrease  with  depth  at  a  rate  defined  by  the  macrofauna 
rework  depth  L\.  For  simplicity,  it  is  assumed  that  the  macrofaunal  and  meiofaunal 
rework  depths  are  equal  ( L\  =  L),  and  the  rate  parameters  for  the  positions  of  the 
source/sink  functions  are  assumed  to  be  identical  (Aa  =  A &). 

The  diffusion  coefficient  and  the  distribution  of  source/sink  events  are  assumed 
to  be  uniform  in  the  horizontal  directions.  In  general,  this  assumption  does  not  seem 
realistic.  However,  it  may  be  reasonable  if  enough  horizontal  averaging  is  performed 
in  the  measurements. 

The  shape  functions  to  be  used  will  be  random  spherical  source  and  sink  func¬ 
tions,  having  the  Fourier  transform 

=  G(k)  =  \  (sinafc  —  ak  cos  ak)  .  (6-14) 

The  radius  a  is  assumed  to  be  a  Gaussian  random  variable  with  normal  distribution 
N(rja;  a a).  Using  equations  (6.12-6.14)  and  the  above  assumptions,  the  spectrum  of 
the  locally  stationary  and  time-dependent  biodiffusion  process  (6.6)  is  rewritten  as 


S„(k, r, r)  = 


(6.15) 


The  expected  value  integral  in  (6.15)  must  be  performed  numerically.  Figure  6.1 
illustrates  an  exponentially  decreasing  distribution  of  spherical  source/sink  functions 
with  random  radii. 
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Top  View 


Horizontal 


Figure  6.1:  Spherical  shape  functions  that  are:  (a)  homogeneously  distributed  in 
the  horizontal;  and  (b)  nonhomogeneously  distributed  in  depth.  The  Poisson  point 
density  in  space  and  time  decreases  exponentially  with  depth  with  a  scale  defined 
by  the  rework  depth  L,  concentrating  activity  near  the  sediment-water  interface. 


The  model  now  consists  of  the  normalized  correlation  coefficient  (6.5)  as  the 
observable  and  the  biodiffusive  spectrum  (6.15)  as  the  forward  model.  The  object 
is  to  minimize  the  difference  between  the  model  predicted  backscatter  correlation 
and  the  data  as  a  function  of  grazing  angle  and  lag  time.  By  using  the  normalized 
correlation,  the  magnitude  of  (6.15)  is  removed  from  the  model  as  an  unknown.  The 
inverse  problem  is  reduced  to  the  solution  of  four  unknowns:  1)  the  mean  radius  r/a 
of  the  spherical  source/sink  shape  functions;  2)  the  variance  cr%  of  the  radius;  3)  the 
biodiffusion  coefficient  Do ;  and  4)  the  rework  depth  L. 

6.3  Orcas  Island  Experiment 

Model  predictions  are  compared  with  data  collected  during  an  acoustic  experiment 
performed  in  a  shallow  water  bay  of  Orcas  Island  in  Puget  Sound  in  Washington 
state.  Acoustic  backscatter  was  recorded  over  an  area  of  seafloor  by  a  bottom 
mounted  transducer  configured  as  illustrated  in  Figure  6.2.  The  acoustic  system  has 
been  described  in  detail  elsewhere  [20,  41],  and  only  a  brief  description  is  given  here. 
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A  calibrated  40  kHz  monostatic  transducer  with  a  fan-shaped  beam  pattern  was 
used  to  scan  a  circular  site  approximately  50  meters  in  radius  at  azimuthal  incre¬ 
ments  of  5°.  A  single  pulse  of  duration  T0  =  2  ms  was  transmitted  for  each  azimuthal 
direction,  and  backscatter  was  recorded  to  produce  a  single  radial  scan.  The  pulse 
has  a  rectangular  envelope,  and  it  was  frequency  modulated  to  cover  a  bandwidth 
of  2  kHz.  The  transducer  was  then  rotated  to  the  next  azimuthal  increment  and 
another  pulse  transmitted.  The  transmission  and  recording  of  backscatter  is  re¬ 
peated  to  complete  a  circular  scan  in  approximately  7  minutes.  A  complete  scan 
was  recorded  at  2.4  hour  intervals  through  the  duration  of  the  experiment.  The 
instrument  remained  at  a  fixed  location  for  approximately  60  days. 

A  total  of  585  circular  scans  was  recorded  to  provide  a  series  of  backscatter  images 
that  reveal  temporal  and  spatial  variability  in  scattering  from  the  sediment.  Figure 
6.3  illustrates  the  spatial  variability  of  backscatter.  The  circular  images  show  the 
mean  and  variance  of  backscatter  estimated  from  50  complete  scans  (over  5  days). 
The  mean  shows  large  spatial  variation  in  backscatter  strength  (>  3  dB).  The  small 
patches  of  high  backscatter  in  the  fourth  quadrant  of  the  circle  are  due  to  targets  and 
manipulations  made  during  the  experiment.  The  second  and  third  quadrants  were 
undisturbed.  The  small  scale  spatial  variations  in  backscatter  (~l-2  meter  scale) 
are  likely  due  to  the  expected  statistical  variations  in  the  signal  (of  ~5  dB).  The  low 
variance  of  the  backscatter  strength  (~  —80  dB)  suggests  that  the  variability  in  the 
mean  backscatter  strength  is  due  to  the  sediment  properties,  rather  than  caused  by 
transients  such  as  reverberation  from  organisms  in  the  water  or  uncorrelated  noise. 
For  both  plots  in  Figure  6.3,  the  Lambert  law  angular  dependence  of  backscattering 
(sin#2)  is  removed.  Both  the  mean  and  the  variance  of  the  scattering  strength  were 
calculated  by  forming  the  ensemble  average  of  the  scattering  cross-section  and  then 
converting  to  dB  using  101oglo. 

The  sediment  in  the  area  is  a  silty-clay  with  moderate  biologically  activity.  No 
bubbles  were  observed  in  the  sediment,  and  an  extensive  set  of  sediment  physi- 


143 


Z 


Figure  6.2:  Scattering  geometry  of  Orcas  experiment.  The  shaded  area  is  a  radial 
scan  for  a  single  pulse  transmission. 


Figure  6.3:  Spatial  variability  in  (a)  the  mean  backscatter  strength,  and  (b)  the 
variance  of  the  backscatter  strength  at  the  Orcas  site.  Data  estimated  from  50 
consecutive  scans  using  40  kHz  backscatter,  and  plotted  on  a  uniform  grid  with  1.5 
meter  spacing. 
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cal  properties  was  measured  at  the  site  by  investigators  of  the  US  Naval  Research 
Laboratory  (NRL).  This  work  included  core  and  in-situ  measurement  of  geoacous¬ 
tic  parameters,  X-radiography  of  sediment  inhomogeneity,  and  stereo-photographic 
measurement  of  sea  floor  roughness.  The  biodiffusion  of  natural  radio-isotope  trac¬ 
ers  was  also  analyzed,  and  an  estimate  of  the  vertical  biodiffusion  coefficient  was 
obtained  in  the  same  location  as  the  acoustic  data  [82]. 

6.3.1  Sediment  Physical  Properties 

The  spectra  of  core  density  fluctuations  were  estimated  using  a  first-order  auto¬ 
regressive  (AR)  model  (see  Appendix  C).  This  model  assumes  that  the  auto¬ 
correlation  of  the  normalized  fluctuations  can  be  expressed  as 

Cpp(r)  =  ^g-inoM/i*  ,  (6.16) 

where  a  is  the  first-order  AR  coefficient  estimated  from  the  core  data  using  the 
Levinson-Durbin  algorithm  [59].  Since  7P  is  a  zero-mean  process,  the  zero-lag  auto¬ 
correlation  is  its  variance,  a2.  The  isotropic  power  law  PSD  function  for  density  fluc¬ 
tuations  is  found  from  (6.16)  where  the  correlation  length  is  defined  as  lc  =  —Sz/  In  a, 
and  Sz  is  the  core  sampling  interval.  The  spectral  parameters  for  the  Orcas  cores 
are  listed  in  Table  6.1  along  with  other  geoacoustic  parameters  supplied  by  NRL.  A 
core  sampling  interval  of  2  cm  was  used  providing  marginal  resolution  compared  to 
the  acoustic  wavelength.  The  core  data  are  also  one-dimensional,  which  necessitated 
an  assumption  of  isotropy  in  the  spectral  estimates. 

6.3.2  Sediment  Volume  Scattering 

Backscatter  strength  per  unit  area  of  sea  floor  was  measured  as  a  function  of  graz¬ 
ing  angle  for  the  Orcas  site.  Owing  to  the  similarity  of  the  magnitude  and  angular 
dependence  of  the  measured  values  compared  to  those  of  other  sites  having  a  range 
of  different  sediment  types  [33,  34,  48],  the  dominant  scattering  mechanism  cannot 
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Density  Ratio 

p/po 

1.406 

Sound  Velocity  Ratio 

c/cQ 

0.977 

Proportionality  Constant 

p 

-0.934 

Compressional  Loss 

8 

0.0019 

Density  Variance 

al 

0.002 

Correlation  Length 

lc 

3.198  (cm) 

Table  6.1:  Sediment  parameters  estimated  from  cores  and  model-data  comparison 
for  the  Orcas  site  [14]. 


be  inferred  from  the  backscatter  alone.  However,  by  using  the  sediment  core  mea¬ 
surements  as  inputs  to  scattering  models,  the  model  results  can  be  compared  with 
observed  backscatter  data,  and  the  dominant  scattering  mechanism  inferred. 

Two  scattering  models  are  considered:  scattering  from  volume  heterogeneity  and 
scattering  from  a  rough  sediment-water  interface.  Surface  roughness  scattering  is 
treated  using  a  first-order  perturbation  approximation  [34]  with  roughness  charac¬ 
terized  by  the  two-dimensional,  isotropic,  power-law  spectrum  estimated  from  sites 
with  similar  sediments  [33,  34].  As  Figure  6.4  shows,  the  predicted  interface  scatter¬ 
ing  was  sufficiently  below  the  measured  values  that  interface  scattering  effects  can 
be  neglected. 

Figure  6.4  shows  the  predicted  backscatter  strengths  found  using  the  first-order 
volume  model  with  and  without  half-space  effects.  In  this  case,  a  first-order  approx¬ 
imation  models  the  data  reasonably  well,  and  half-space  effects  are  negligible  over 
the  angular  range  of  the  data.  The  difference  between  model  predictions  and  data 
is  most  likely  due  to  errors  in  the  estimated  heterogeneity  spectra. 

No  significant  oceanographic  processes  or  sediment  transport  events  were  ob¬ 
served  during  the  experiment.  For  the  Orcas  site,  sound  speed  fluctuations  were 
recorded  for  the  duration  of  the  experiment  by  conductivity,  temperature,  and  depth 
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Figure  6.4:  Comparison  of  measured  backscatter  with  model  predictions  for  the 
Orcas  site. 


i 
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Figure  6.5:  Sound  speed  fluctuations  in  the  water  during  the  Orcas  experiment. 


(CTD)  measurements  (see  Figure  6.5).  Given  the  nature  of  the  sediment  and  the 
expected  low  level  of  roughness  scattering,  it  is  expected  that  acoustic  penetration 
into  the  sediment  and  the  effects  of  biological  mixing  are  the  dominant  mechanisms 
for  decorrelation. 

6.3.3  Model-Data  Comparison 

The  model  for  the  temporal  correlation  of  seafloor  backscattering  due  to  biodiffu-  , 
sive  sediment  mixing  was  compared  with  backscatter  correlation  estimated  from  the 
Orcas  experiment..  Figure  6.6  shows  a  comparison  of  the  observed  and  predicted 
correlations  as  functions  of  lag  time.  The  model  was  fitted  to  the  data  (plotted  as 
points)  by  trial  and  error  while  constraining  the  parameters  of  the  sediment  model 
to  realistic  values.  The  offset  of  the  observed  correlation  profile  from  the  predicted 
profile  is  due  to  uncorrelated  noise  in  the  data.  Reverberation  in  the  water  column 
that  is  uncorrelated  from  scan  to  scan  is  a  likely  source  of  noise. 

Table  6.2  outlines  the  parameters  estimated  from  the  model-data  comparison. 
The  estimated  biological  parameters  are  typical  for  the  type  of  sediment  and  envi¬ 
ronment  encountered  at  the  Orcas  site.  An  estimate  of  the  biodiffusion  coefficient 
for  naturally  occurring  radioisotopes  (A)  was  measured  from  cores  in  the  same 
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Diffusion  Coefficient 

Do 

2.4  [cm2 /year) 

Macrofauna  Rework  Depth 

Lx 

3.8  (cm) 

Meiofauna  Rework  Depth 

L 

3.8  (cm) 

Mean  Radius 

*?a 

0.4  (cm) 

Radius  Standard  Deviation 

cra 

0.14  (cm) 

Table  6.2:  Sediment  parameters  estimated  from  cores  and  model-data  comparison 
for  the  Orcas  site. 


area  as  the  acoustic  experiment.  The  core  analysis  produced  an  estimate  of  Db  «  8 
cm2/year  [82].  This  value  is  higher  than  the  acoustically  inferred  value  of  Do  ~  2.4 
cm2 /year.  One  explanation  is  that  the  value  estimated  from  the  cores  is  a  different 
type  of  diffusion  coefficient  than  the  coefficient  estimated  using  (6.12).  The  D  found 
using  the  bioturbation  model  developed  in  this  thesis  is  the  meiofaunal  diffusion 
coefficient.  The  Db  estimated  from  cores  uses  the  classical  diagenesis  model  of  dif¬ 
fusive  mixing,  and  therefore  it  includes  both  mieofaunal  and  macrofaunal  mixing. 
In  general,  it  would  be  expected  that  the  meiofaunal  diffusion  coefficient  is  smaller 
than  the  total  diffusion  coefficient. 

6.4  Summary  and  Conclusions 

A  biodiffusion  model  using  random  excitation  to  describe  the  time  evolution  of 
sediment  inhomogeneities  is  combined  with  the  first-order  perturbation  treatment 
of  sediment  volume  backscatter.  The  first-order  perturbation  model  for  volume 
backscatter  is  validated  for  use  at  the  Orcas  site  by  data-model  comparison.  Half¬ 
space  effects  are  shown  to  be  unimportant  at  the  low  grazing  angles  used  in  the  Orcas 
experiment.  Sediment  physical  properties  and  inhomogeneity  spectra  were  estimated 
using  measured  core  data.  Higher  resolution  two-  or  three-dimensional  core  analyses 
(e.g.,  X-radiographs,  CT)  are  recommended  for  future  work  to  improve  the  sediment 
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Figure  6.6:  Comparison  of  measured  backscatter  correlation  and  model  predictions, 
(a)  9  =  14.7°,  (b)  9  =  9.4°,  (c)  0  =  7.6°,  (d)  9  =  6° 
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ground  truth  estimates. 

Backscatter  measured  at  the  Orcas  site  was  used  to  estimate  the  temporal  cor¬ 
relation  of  the  scattered  field.  The  correlation  estimate  was  then  used  to  estimate 
benthic  biological  activity  by  comparison  with  the  biodiffusive  mixing  model.  Us¬ 
ing  the  correlation  data  for  the  Orcas  site,  the  proposed  model  of  sediment  mixing 
predicts  sediment  biodiffusion  parameters  that  are  realistic.  Verifying  the  acousti¬ 
cally  inferred  biodiffusive  parameters  is  problematic  because  there  is  no  analogous 
sampling  method  that  provides  similar  spatial  and  temporal  resolution.  This  is  the 
main  reason  acoustic  remote  sensing  methods  are  interesting.  Further  model-data 
comparisons  and  better  ground  truth  estimates  of  sediment  physical  and  biological 
properties  are  needed  to  verify  the  methodology.  However,  this  initial  and  simple 
application  of  the  model  to  data  provides  encouraging  results. 
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Appendix  A 

HALF-SPACE  GREEN’S  FUNCTION 

Consider  a  two-dimensional  medium  with  a  source  at  position  r'  =  ( x ',  z')  and 
a  receiver  at  position  r  =  (x,  z)  both  in  the  LHS  of  two  homogeneous  half-spaces. 
The  Green’s  function  is  the  solution  of  the  wave  equation, 

[V2  -|-  P]  #)(r,  r')  =  -5(r  -  r')  ,  (A.l) 

that  satisfies  reciprocity  and  the  appropriate  boundary  conditions  along  the  interface 
(at  z  —  0).  The  solution  is  found  by  first  taking  the  Fourier  transform  in  x  to  obtain 
the  spectral  representation  of  (A.l), 

+  Sb  =  -  *0  ,  (A-2) 

where 

t>  =  s/K-kl  ■  (A.3) 

The  solution  for  (A.2)  in  the  region  z  <  0  and  z  >  z'  is 

g0  =  A  [R{kx)e~ipz  +  eipz]  , 
and  in  the  region  z  <0  and  z  <  z'  is 

go  =  Be-**’  . 

The  coefficient  R  is  the  reflection  coefficient  between  the  two  regions: 

R  __  PoP_  —  0o 
Po(3  +  pfi 0 


(A.4) 
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The  unknown  coefficients  A  and  B  are  found  by  matching  the  solutions  and  imposing 
the  discontinuity  in  the  normal  derivative  at  z  =  z'  required  by  (A. 2),  such  that 


A  ^R(kx)e~ipz'  +  eipz'j  =  Be~ipz' 


and 


if3A  \—R{kx)e~ipz>  +  eipz'\  +  if3Be~ipz'  =  -^-e~ikxX'  . 
L  J  2tt 

Solving  for  A  and  B, 

A  —  l  -ikxx'  -ij3z' 

A  47 r/3 


and 


5  =  4 '-f***' 


e-ipz'  +  eipz> 


the  combined  solution  for  the  entire  region  £  <  0  is 


go  = 


4k  (3 


e-ikxx'  ^R(-k^e-il3(z+z>)  +  ei0\z-z' |j 


The  final  solution  is  found  by  taking  an  inverse  Fourier  transform  in  x: 

<7o(r;r')  =  ^  [R(kx)eip^  +  eip^  . 

The  Green’s  function  for  a  receiver  in  the  UHS  (z  >  0)  and  a  source  in  the  LHS 
{z!  <  0)  is  found  in  a  similar  manner.  The  solution  for  (A. 2)  in  the  region  z  >  0  is 

go  =  Ceip°z  , 


where 

/So  =  \Ao2  -  ki  •  (A-S) 

Again,  solve  for  the  unknown  coefficients  C  by  matching  the  solutions  and  their 
derivatives  at  z  =  0,  such  that 

C  =  A(R  +  1)  =  AT’ , 
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and 

Hkc  =  'lA(i  -  r)  . 

Po  P 

The  coefficient  T'  is  the  transmission  coefficient  from  the  LHS  to  the  UHS: 

T,  _  2PoP 

poP  +  pPo 

Using  the  derivation  of  A  from  the  reflection  Green’s  function  above,  the  solution 
for  the  spectral  component  of  the  transmission  Green’s  function  is 


9o  = 


4wp 


T'(kx)e 


ikxx *  gi/3 o  z—ij3zf 


The  final  solution  is  found  by  performing  an  inverse  Fourier  transform  is 
<?o(r;r')  =  ±  J°°  ^.T' {kx)eik^x-x^ . 
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Appendix  B 

FAR-FIELD  HALF-SPACE  GREEN’S  FUNCTION 


The  half-space  Green’s  function  is  evaluated  in  the  far-field  of  the  interface  using 
the  method  of  stationary  phase  or  saddle-point  method  [29].  Begin  with  the  two- 
dimensional  Green’s  function  for  transmission 


OO 

g0( r,r')  =  *  J  T'^x)  e^'^'  eik^x~x'Ukx  , 

— OO 

(B.l) 

where 

P=yJW-kl, 

/ - 

(B.2) 

ii 

<SE 

1 

and 

Tf  ZpoP 

(B.3) 

pofl  +  PP o 

Change  integration  variable  in  (B.l)  to  cylindrical  coordinates  with  kx  =  ko  cos  9  and 
f30  =  k0  sin  6,  and  define  the  the  position  r  =  (x,  z)  as  x  =  r  cos  9S  and  z  =  r  sin#s. 
Then  using  the  identity 


cos(0  -  9S)  =  cos(0)  cos(0s)  +  sin (9)  sin(0s)  ,  (B.4) 


the  integral  (B.l)  can  be  rewritten  as 

»(r,0  =  ^jfr(d)^ 


(B.5) 


The  complex  contour  C  is  defined  in  Figure  B.l.  The  function  F{9)  is  defined  as 

(B.6) 


F(0j  _  _r'(^)e-it„'eo,<-,t„'.ln»|toSin^  1 
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Re[0] 


Figure  B.l:  Complex  contour  in  the  6  plane. 


and  the  exponent  function  f{6)  is  defined  as 


f(9)  =  ik0cos(8  —  6S)  .  (B.7) 

It  can  be  shown  that  the  integrand  vanishes  at  both  ends  of  the  contour  [29]. 

The  stationary  phase  point  is  given  by 

dm 


86 


=  —  ikosm(6  —  6S )  =  0, 


(B.8) 


or  6  =  0S.  The  function  f{6)  is  then  expanded  about  6S  as 

(« -  «.)2  „„ 


m = /(«.) + 


■/"(«.)  , 


(B.9) 


where  the  double  primes  denote  the  second  derivative  with  respect  to  6.  The  in¬ 
tegrand  can  then  be  approximated  as  a  Gaussian  curve,  which  for  large  r  can  be 
evaluated.  Evaluating  the  integral  gives  an  approximate  expression  for  the  far-field, 

1 1/2 


<7o(r,r')  =  — 


27T 


(B.10) 


Y\S"(».)\r\  ’ 

where  7  =  — 7t/4  is  chosen  to  properly  orient  the  path  of  integration  in  the  complex 
plane. 
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The  far-held  Green’s  function  is  then  written  as 

nJv  T'(0S)  tfc0r— ak,-r'  — t7r/4 

ffo(r’r,_  f)  s/W^  ’ 

or  equivalently  as 

nJr-  r'\  —  i—  pik0r  -iks-r'  -i-ir/4 

9“(r’r,“  VvC5?e  6  e  ’ 

where 

T_  tyPo 
PoP  +  ppo 

The  same  procedure  is  performed  to  find  the  three-dimensional  far-field 
function: 


0o(r,r') 


PoTtysY  ikor  —iks'r* 
p4nr 


(B.ll) 

(B.12) 

(B.13) 

Green’s 


(B.14) 


Appendix  C 


AUTO-REGRESSIVE  SPECTRAL  ESTIMATION 


Auto-regressive  spectral  estimation  is  based  on  modeling  the  digital  data  se¬ 
quence,  y(n),  as  the  output  of  a  linear  system, 


y(n )  =  h(n)  *  x(n )  , 


(C.l) 


characterized  by  the  rational  system  function 


H(z) 


1 

W) 


i  +  Yh<tkz  k 

k=l 


(C.2) 


If  the  input  and  output  of  the  system  are  stationary  random  processes,  the  auto¬ 
correlation  of  the  data  sequence  is  defined  in  terms  of  the  input  auto-correlation  and 
the  discrete  step  response  of  the  system: 


Ryy(m)  =  Rxx(m)  *  h(m)  *  h*(—m )  .  (C.3) 

Using  an  estimate  of  Ryy(m)  taken  from  the  data  and  assuming  the  input  is  a  zero- 
mean  white  noise  process  with  auto-correlation  Rxx(m)  =  q8(m ),  determining  the 
spectrum  of  the  output  process  becomes  a  problem  of  estimating  the  coefficients  a*,. 
The  constant  q  is  the  variance  of  the  input  process.  Several  well  known  techniques 
such  as  the  Yule- Walker  method  and  the  Levinson-Durbin  algorithm  can  be  used 
to  estimate  these  parameters  [59]  [55].  With  an  estimate  of  the  system  function  (or 
system  filter  coefficients),  the  spectrum  of  the  output  process  is  found  by  evaluating 
equation  (C.2)  on  the  unit  circle, 
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where  z  —  etk  and  k  is  the  spatial  frequency. 

A  first-order  AR  process  is  characterized  by  a  system  function  with  a  single  pole 

H(z)  =  (C.5) 

1  —  az  1 

For  a  <  1  the  impulse  response  function  of  this  system  is 

h(n)  =  anu(n )  .  (C.6) 

The  auto-correlation  of  this  process  is  found  using  equations  (C.3)  and  (C.6)  to  be 

*Um)  =  (°-7) 

The  normalized  auto-correlation  function  is 

ryy(m )  =  a|m|  =  (C.8) 

In  general,  we  want  to  consider  the  correlation  function  of  a  continuous  process  in 
3-D  space.  Thus,  if  we  generalize  equation  (C.7)  to  represent  an  isotropic  process 
in  which  m  =  r  =  -sjx2  +  y2  +  z2  and  apply  the  Wiener-Khinchin  theorem  [30],  the 
spatial  spectrum  of  a  first-order  AR  process  is  found  to  be 

Syy(k )  =  (a2+^2)2  ’  (C-9^ 

where  a  =  —  In  a/ 5m,  (3  =  q/(  1  —  a2),  and  8m  is  the  sampling  interval. 

How  does  this  relate  to  power  law  spectrum  of  normalized  sediment  density 
fluctuations  derived  from  an  exponential  correlation  function  (as  in  Section  2.1.1)? 
By  the  definition  of  (C.8),  a  =  l//c.  The  variance  of  the  normalized  sediment 
density  fluctuations  is  simply  the  variance  of  the  input  process,  a2  =  q/(  1  —  d2). 
Thus,  equation  (C.9)  can  be  written  as 

<7^ 

S”W  =  n%(l /£+k*y  ' 


(C.10) 
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Appendix  D 

FFT  EVALUATION  OF  THE  HALF-SPACE  GREEN’S 

FUNCTION 

Since  the  unknown  field  points  are  discretized  on  a  uniform  grid,  the  wavenumber 
integration  of  the  half-space  Green’s  function  (2.33)  can  be  performed  efficiently 
using  an  FFT  [38].  First,  use  the  discretization  of  the  x-coordinate, 

x„  —  •Emin  +  fiAs  ,  for  Tt  —  1, 2,  •  ■  .  ,  N  ,  (HT ) 

to  define  the  corresponding  discretized  wavenumber  in  the  x-direction, 

kXm  =  "f  ,  for  TTi  —  1,  2,  •  •  •  ,  ilf  ,  (D.2) 

where 

AzA k,  =  -  .  '  (D.3) 

Substituting  xn  and  kxm  into  the  definition  of  the  half-space  Green’s  function  for 
reflection  and  approximating  the  wavenumber  integration  as  a  summation, 

J  dkx  (D.4) 

the  Green’s  function  is  rewritten  as 

°°  AT, 

^  iS  P»  (D.5) 

\zn  ~ Z*  |  j  ^,rn (&min 

With  attenuation  in  the  medium,  the  term  in  square  brackets  will  vanish  at  high 
wavenumbers  and  the  summation  can  be  truncated.  The  summation  now  takes  the 
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form  of  an  inverse  FFT, 

{R(hxmyl3m\zn+z>\  _j_  eipm\zn-z'\^eim&kx(xmin-x')  , 

fim 

(D.6) 


iM 

9o(xm  zn]  x',  z')  =  —FFT 


-1 


where 

1  M- 1 

FFT-1  [F(x„,  2„)]  =  F(xm2„)e"'T  . 

771=0 


(D.7) 


The  computation  of  the  Green’s  function  is  greatly  improved.  A  single  FFT 
computation  is  required  to  find  a  row  of  values.  The  discretization  spacing  in  the 
^-direction  must  be  chosen  to  provide  adequate  spacing  in  the  spectral  domain  to 
fully  represent  the  spatial  frequency  characteristics  of  the  integrand.  If  Ax  or  M  is 
too  small,  aliasing  may  corrupt  the  FFT.  Owing  to  the  assumed  periodicity  of  the 
FFT,  the  range  over  which  the  function  F  is  evaluated  must  be  sufficiently  large  so 
that  F  approximately  vanishes  outside  the  interval. 
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Appendix  E 


GENERATING  GAUSSIAN  RANDOM  MEDIA 


Realizations  of  a  two-dimensional  Gaussian  random  medium  are  generated  by 
discrete  inverse  Fourier  transforms.  Consider  the  random  variable  f(x,  z)  and  the 


discrete  Fourier  transform 


f{x,z)  = 


EE* 


ikXnx 

77171 c  c 


(E.l) 


where  kzm  =  2irml  Lz  and  kxn  =  2im/Lx.  The  correlation  of  f(x,  z)  is  defined  as 


Rff(x  1  -  x2,z !  -  z2 )  =  (f(xuzi)f*(x2,z2)) 


(E.2) 


If  f(x,  z)  is  a  Gaussian  random  variable,  then  the  coefficients  bmn  are  independent 
and  uncorrelated,  and  they  can  be  defined  as  the  Gaussian  random  variable 


6mn  =  ^[iV(0,1)+*W(0,1)1  ’ 


(E.3) 


where  Amn  is  the  amplitude.  The  function  N( 0, 1)  is  a  uncorrelated  Gaussian  random 
number  with  zero  mean  and  unit  variance,  such  as  a  numerical  random  number 
generator.  Using  (E.l)  the  correlation  of  f(x,z)  can  be  rewritten  as 

D  /~  „  ~  _  \  _  (^7I’)  \'\  '  D  eikzm(zi-z2) eikxn{xi-x2)  (T?  A  \ 


Rjf{ Xl  ~  X2,Zi  -  Z2)  = 


(LxLzy 


EE* 


where 


Bmn  —  ( |  bmn  |  ) 


(E.5) 


The  correlation  of  f(x,  z )  is  alternatively  expressed  as  the  discrete  Fourier  trans¬ 
form  of  the  spectrum  of  the  process  that  generated  f(x,z ): 
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Equating  (E.6)  and  (E.4),  the  random  Fourier  coefficients  of  the  correlation  function 
are  related  to  the  spectrum  as 

=  (E.7) 

(27rr 

Now,  using  (E.3)  and  (E.5)  in  (E.7),  the  amplitudes  of  the  random  Fourier  coeffi¬ 
cients  of  f(x,  z)  are  found  to  be 


A  -  J- 

mn  (2  tt) 

Using  a  desired  spectrum  for  S '//,  such  as  the  power-law  spectrum  (2.18),  and  gen¬ 
erating  the  Fourier  coefficients  using  (E.3)  and  (E.8),  a  random  realization  of  f(x,  z) 
can  be  generated  by  performing  the  inverse  Fourier  transform.  A  discrete  random 
realization,  f(xi,zj),  is  found  by  discretizing  the  finite  region  (LX,LZ)  as  a  uniform 
grid  and  performing  the  inverse  transform  at  the  grid  locations: 
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Appendix  F 

MASS  OPERATOR  FOR  THE  BILOCAL 
APPROXIMATION 

The  mass  operator  for  Dyson’s  equation  using  the  bilocal  approximation  (4.17) 
is  defined  in  the  wavenumber  domain  as 

M( k)  =  k4MKK(k)  -  2iPk  ■  MKp(k)  -k-Mpp(k)-k,  (F.l) 


where 


Mkk( k)  =  n2  J  Cpp(r)g0(r)e  ik  rdr  , 

(F.2) 

M«p(k)  =  n  J  Cpp(r)Vg0(r)e~tk  rdr  , 

(F.3) 

and 

Mpp(k)  =  J  Cpp(r)VVg0(r)e-ik-rdv  . 

(F.4) 

Using  the  free  space  Green’s  function 

Akr 

*<r>  “  5? 

(F.5) 

with  an  isotropic  and  exponential  correlation  function 

C„( r)  =  <rp2e-W/'-  , 

(F.6) 

each  mass  operator  can  be  derived  analytically. 

The  integrals  of  (F.2-F.4)  are  solved  by  changing  integration  variables  to  spherical 
coordinates  and  orienting  k  in  the  ^-direction,  so  that  k  •  r  =  kru  and  u  —  cos  9. 
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Since  there  is  only  one  wave  vector  k  in  (F.l),  and  it  is  pointed  in  the  ^-direction, 
only  the  z-component  of  the  vector  and  the  ^-component  of  the  dyad  ML,/;  are 
considered.  Therefore,  equations  (F.2)  and  (F.3)  are  rewritten  as 


//2/t2  roo  /*1 

MKK(k)  =  ^P  /  e-{l/u-ik+iku)rrdrdu 

2  Jo  J- i 

(F.7) 

and 

k  •  MKp(k)  =  ^  J°°  J'  [ikr  -  1]  e-Wl‘-ik+ikuKdudr  , 

(F.8) 

using 

d  dr  d  d 

dz  dzdr  U  dr 

(F.9) 

Through  integration  by  parts,  the  dyadic  mass  operator  (F.4)  is  rewritten  as 

r  1  OO 

k •  Mpp(k)  •  k  =  |fc2Cpp(r)e“,kr^o(r) 

-  k2  J  j-z  lcpp(r)e~ik'r }  §;So(r)  dr  , 

(F.10) 

which  after  applying  the  radiation  condition  and  spherical  coordinates,  is  rewritten 
as 


L2-2  roo  ri  _ 

k  •  M„(k)  •  k  =  — /  /  \ikr  -  1]  [«  +  ikle)  e-Wl‘-ik+iku>ududr  .  (F.ll) 

2 lc  Jo  J-  i 

The  integration  on  r  and  u  in  all  three  integrals  can  now  be  performed  to  obtain 

a-  2/2 

(F.12) 


(F.13) 


(F.14) 
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where 


a  =  1  —  iklc  , 
b  =  iklc  , 


(F.15) 


c  =  ln[(a  +  b)/(a  —  6)]  . 
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